© The Authors 2026

1. Introduction

Millisecond pulsars (MSPs) represent the oldest objects in the entire population of pulsars. Although they are at the very end of stellar evolution, they are the fastest spinning pulsars known (with periods P < 30 ms). In contrast to normal isolated pulsars and magnetars, MSPs are found in binary systems, a fact that provides strong evidence of the rapid rotation. During the binary evolution stage, they likely accreted matter from their companion, which increased their angular momentum and spun them up to high rotational frequencies (Alpar et al. 1982; Bhattacharya & van den Heuvel 1991; Tauris et al. 2012). More broadly, these objects are referred to as ‘recycled pulsars’, as they are thought to have initially undergone a significant spin-down before being spun up again. During the accretion phase, these sources have been observed as X-rays binaries (Wijnands & van der Klis 1998; Patruno & Watts 2021). Although different mechanisms (accretion process, magnetic field decay) are involved in recycling of pulsars, they are not completely understood. Nevertheless the recycling scenario is widely accepted in the community today (Radhakrishnan & Srinivasan 1982; Yang & Li 2023; Tauris & van den Heuvel 2023). MSPs without a companion have also been observed, and their properties compared to recycled pulsars with a companion are not distinguishable, which supports the idea that all recycled pulsars must have been related to binary evolution at some point (Srinivasan & van den Heuvel 1982). Moreover, recycled pulsars appear in various forms. Some are isolated, as mentioned earlier (Backer et al. 1982); others have a companion that can be degenerate or not; and in some cases the MSP in the binary can be considered a spider. These include black widows, where the pulsar’s very light companion (M < 0.05 M_⊙) is being destroyed by the pulsar wind, and redbacks, which have more massive companions (M ∼ 0.1 − 0.4 M_⊙) that are ablated by the pulsar (Roberts 2013; Blanchard et al. 2025). Furthermore, transitional MSPs are sources that switch between X-ray and radio emitters (Papitto & de Martino 2022). In eclipsing MSPs, the radio emission is eclipsed by the companion and could be due to ejected ionised material from the companion that is absorbing the radio waves from the MSP (Fruchter et al. 1988). Finally, MSPs with planets also exist (Wolszczan & Frail 1992).

Until 2008, very few pulsars had been detected in γ-rays (about seven; Thompson 2008), and it is thanks to the Fermi Gamma-Ray Space Telescope (Fermi) that this number was significantly increased to about 300 today (Smith et al. 2023). Only a few MSPs have been detected in gamma-rays without a radio counterpart. Already in the 1990s, it was proposed that many MSPs might be γ-ray sources (Srinivasan 1990), and Fermi has since led to a dramatic increase in the number of detected γ-ray MSPs, although the vast majority of known MSPs have been discovered through radio observations. Radio pulsar surveys, such as those conducted with the Parkes, Arecibo, and Green Bank telescopes; the Five-hundred-meter Aperture Spherical Telescope (FAST); and MeerKAT, have identified over 400 MSPs to date (Manchester et al. 2005; Qian et al. 2020; Bailes et al. 2020), many of which were later found to emit in γ-rays as well. Radio observations remain essential for accurately determining pulsar timing parameters, detecting binary companions, and characterising their orbital dynamics. Furthermore, several MSPs are still only visible in the radio domain, highlighting the complementarity of multiwavelength studies.

A powerful tool for investigating neutron star (NS) evolution is population synthesis. Numerous studies have applied this method to different NS families, including the normal pulsar population (Faucher-Giguère & Kaspi 2006; Gullón et al. 2014; Sautron et al. 2024; Pardo-Araujo et al. 2025), magnetars (Jawor & Tauris 2022; Beniamini et al. 2019; Sautron et al. 2025), and MSPs (Kiel et al. 2008; Story et al. 2007; Gonthier et al. 2018; Tabassum & Lorimer 2025). From the simulations in these studies, invaluable insights can be extracted regarding the underlying population properties such as the birth spin period, magnetic field decay timescale, and birth rate as well as the emission processes at play. In this approach, NSs are generated at birth and evolved to the present epoch, after which detection criteria are applied to determine whether they would be observable with current or future instruments.

In this work, we focus exclusively on mildly or fully recycled pulsars (a mildly recycled pulsar has accreted matter but not enough to become fully recycled with P < 30 ms, and thus it has 30 ≤ P ≤ 1000 ms and a Ṗ < 10⁻¹⁶) that are not in globular clusters (GCs). Although about half of the known population of MSPs belong to a GC, studying pulsars that are in GCs is very different. Indeed, to study them accurately with simulations, we would need to take many interactions into account. Since there are high stellar densities in these environments (much higher than the Galactic field), we would need to take into account the fact that there could be an exchange of companions (at least it should happen more compared to the Galactic field), there could be triple systems, or there might be a different distribution at birth for the eccentricity or the semi-major axis of the binary. For instance, Berteaud et al. (2024) used simulation-based inference (Cranmer et al. 2020) to study the MSP population in GCs, and they focused on Terzan 5 by finding constraints on the luminosity function of its MSPs and the number (with uncertainties) of MSPs that Terzan 5 hosts.

Story et al. (2007) and Gonthier et al. (2018) are the most complete studies about the population synthesis of MSPs. The less recent study of Story et al. (2007) made predictions for Fermi and revealed that many MSPs would be detected by it as radio quiet γ-ray sources, which was indeed true (Abdo et al. 2013; Smith et al. 2019, 2023). In their paper, the evolution of the MSPs started after their spun-up episode was finished, the magnetic field of MSPs was constant, and only MSPs from the Galactic disc were considered. The more recent study of Gonthier et al. (2018) improved on the studied sample by considering MSPs in the bulge, too. They predicted that 11 000 MSPs would be necessary in the Galactic bulge to explain the γ-ray galactic center excess detected in the Galactic centre, and their results are in agreement with those derived from Fermi detections. Although the results of Gonthier et al. (2018) are compelling, our approach, which self-consistently models the full binary evolution alongside the pulsar’s evolution and a detection model in the γ-ray and radio to reproduce the recycled population, offers several advantages. Our study is so far the only one that allows to follow the full evolution to the MSP state and even to produce mildly recycled pulsars while taking into account the detectability in both radio and γ-ray wavelengths. Our study constitutes a valuable opportunity to investigate whether the accretion-driven formation channel is an effective pathway for producing MSPs, especially in comparison to alternative scenarios such as the collapse of white dwarfs (WDs) into rapidly rotating NSs. Moreover, it enables the derivation of statistical properties of the recycled population, which can be otherwise difficult to access observationally. Finally, it allowed us to make predictions for upcoming surveys, such as the Square Kilometre Array (SKA¹).

This paper explores many aspects of the binary evolution leading to recycled pulsars. It is organised as follows. In Sect. 2 we describe the model to generate and evolve the pulsars. We describe the model to detect the pulsars in the simulation in Sect. 3. Finally, the results are presented in Sect. 4, discussed in Sect. 5, and summarised in Sect. 6.

2. Evolution model

Compared to previous studies, the novelty of the present work lies in the careful treatment of the evolution of binary systems up to the formation of compact objects in their final stage. At the beginning of the simulation, all of the secondary stars considered in the binaries are main sequence (MS) stars that are at a certain percentage of their lives. This percentage is taken from a random uniform distribution between 0 and 100%, while the NS in the binary starts at age 0. We begin by presenting the stellar evolution code, Stellar EVolution for N-body (SEVN), followed by a description of the binary processes it includes. We then introduce the adopted birth distributions and initial characteristics, and outline the prescriptions implemented for the evolution of NSs. Next, we discuss the spatial evolution of the systems, before finally addressing the concept of the death valley.

2.1. The SEVN code

The SEVN code is a state-of-the-art binary population synthesis code. It models single stellar evolution by interpolating pre-computed stellar tracks on the fly, and accounts for binary interactions using analytic and semi-analytic prescriptions (Spera & Mapelli 2017; Spera et al. 2019; Mapelli et al. 2020; Sgalletta et al. 2023). In this work, we used the version of SEVN described in Iorio et al. (2023)². In Sects. 2.2, 2.3, 2.4, and 2.5 we detail a general overview of the main features that we used from SEVN, with modifications on the NS evolution proposed by SEVN to get the best description possible. We refer to Iorio et al. (2023) for a detailed description of the code. The SEVN code requires a maximum NS mass that can be reached in the simulation to be set. Thus, we selected this maximum to be 3 M_⊙, as in the work of Sgalletta et al. (2023).

2.2. Binary processes

Various binary processes come into play during the course of binary evolution, in the following we describe the ones relevant in this work that are implemented in SEVN. The Roche lobe R_L, of a star in a binary system defines the region of space within which matter remains gravitationally bound to the star. When a star fills its Roche lobe, material can flow towards the companion star under the influence of its gravitational pull, a process known as Roche-lobe overflow (RLO). The RLO affects several key parameters, including the mass ratio, the individual masses and radii of the stars, and the binary’s semi-major axis. At each time step, SEVN evaluates the Roche-lobe radii of both stars using the analytical formula derived by Eggleton (1983):

$$\begin{array}{c}\hfill \frac{R_L}{a}=\frac{0.49q^{2/3}}{0.6q^{2/3}+ln(1+q^{1/3})}·\end{array}$$(1)

If the radius of the companion of the NS satisfies r ≥ R_L, then an RLO episode starts and mass falls from the secondary star, which filled its Roche lobe, to the NS (the opposite case is also possible, although very unlikely as the system would need to be very compact). The mass transfer can then either be stable or unstable. In order to decide in which case the system is, in SEVN when the mass ratio q is greater than a critical value q_c that depends on the stellar evolutionary phase, the mass transfer is unstable on a dynamical timescale. The values of q_c used are the same as in Hurley et al. (2002), Spera et al. (2019), Giacobbo & Mapelli (2020); see the values in the first column of Table 3 in Iorio et al. (2023).

A stable mass transfer results in the mass-loss rate, Ṁ_d, proposed by Hurley et al. (2002):

$$\begin{array}{c}\hfill \dot{M_d}=-F(M_d){(ln\frac{R_d}{R_{L,d}})}^3{\mathrm{M}}_{\odot }{\mathrm{yr}}^{-1},\end{array}$$(2)

where F(M_d) is a normalisation factor, R_d is the radius of the donor and R_L, d is the Roche lobe of the donor star. SEVN accounts for non-conservative mass transfer, meaning that the donor star can lose more mass than the companion actually accretes. The accreted mass, Ṁ_a, is then modelled according to the following prescription:

$$\begin{array}{c}\hfill \dot{M_a}=\{\begin{array}{c}\min ({\dot{M}}_{\mathrm{Edd}},-f_{\mathrm{MT}}\dot{M_d})\mathrm{if}\mathrm{the}\mathrm{accretor}\mathrm{is}\mathrm{a}\mathrm{compact}\mathrm{object}\hfill \\ -f_{\mathrm{MT}}\dot{M_d}\mathrm{otherwise},\hfill \end{array}\end{array}$$(3)

where Ṁ_Edd is the Eddington accretion rate and f_MT is the mass accretion efficiency, which is within the range [0,1]. In SEVN, f_MT = 0.5 (Bouffanais et al. 2021; Iorio et al. 2023).

In the case of unstable mass transfer, the binary results in either a merger or a common envelope (CE). There are debates in the community about if it is possible to have spin-up during the CE phase (Osłowski et al. 2011; MacLeod & Ramirez-Ruiz 2015; Chamandy et al. 2018; Chattopadhyay et al. 2020). Because of the controversy, spin-up during CE is not included in SEVN. Therefore, if the mass transfer is not stable, in this simulation the pulsar cannot be recycled, and that is why we do not describe the CE evolution here (but see Iorio et al. 2023; Sgalletta et al. 2023 for more details).

2.3. Birth properties

Many birth characteristics generated for each pulsar are independent of their companion, such as the position in the Galaxy (x₀,y₀,z₀), inclination angle (χ₀), spin period (P₀), magnetic field (B₀). Therefore, these quantities are generated from the same distributions, with the same parameters, as described in the paper of Sautron et al. (2024). However, compared to our study of the normal pulsars, we did change two elements in this work: firstly, the kick velocity distribution was modified by adopting a lower dispersion of σ_v = 70 km s⁻¹ which better matches the recycled pulsar population and is consistent with earlier works (Story et al. 2007; Gonthier et al. 2018). Indeed, this prescription is to reflect the assumption that the binary survived the supernova explosion of the primary star that formed the NS. All simulated systems are thus assumed to remain bound after the first supernova event that initially formed the NS in the binary. Furthermore, this prescription provides a better match to the recycled pulsar population, as a high velocity dispersion would result in the loss of a large fraction of systems (resulting in ∼3.5 times fewer detectable NSs), thereby requiring a significantly higher birth rate to reproduce the observed population. Nevertheless, if the secondary star undergoes a supernova explosion to become an NS or a black hole (BH) later in the evolution, the binary may still be disrupted due to the associated mass loss and natal kick from the formed NS or BH, as implemented in the SEVN code. Secondly, instead of using the radial distribution in the galactic plane of Yusifov & Küçük (2004), we used a truncated normal distribution (the distribution is truncated to be always greater than 0) with a mean μ_R = 11.3 kpc, and a standard deviation σ_R = 1.8 kpc to have the NSs born in the spiral arms (see Sautron et al. 2024 to check how the NSs are then placed into the spiral arms). Actually, using the radial distribution proposed by Yusifov & Küçük (2004) causes too many pulsars in the Galactic centre to be detected compared to observations. That is why we opted for the truncated normal distribution for the radial distribution in the spiral arms, which was also an option in Faucher-Giguère & Kaspi (2006) and Tabassum & Lorimer (2025).

The population of NS generated have their ages taken from a random uniform distribution between a minimum age of 5 × 10⁷ yr, as it is very unlikely that a pulsar can be recycled within a shorter time, and a maximum age of 13.9 × 10⁹ yr, which is within the range of the estimated age of the universe (Riess et al. 1998). In each simulation, the number of sources is characterised by a birth spacing X, where X is chosen as an integer number of years between the births of subsequent pulsars. In this paper, we report this quantity as either the birth spacing X in years or the birth rate R_recycled = (1/X) yr⁻¹. This approach can be used to verify whether a constant birth rate is realistic, by comparing the results from the simulation to observation. Before running the simulation, we choose the number of binaries, and then by computing $R_{\mathrm{recycled}}=\frac{\mathrm{Recycled}\mathrm{pulsars}\mathrm{formed}}{{\mathrm{age}}_{\max }-{\mathrm{age}}_{\min }}$ at the end of the simulation, we get the birth rate. In this work, simulating between 4.6 × 10⁵ and 5.7 × 10⁵ pulsars in a binary allowed us to get similar P−Ṗ diagram between observations and simulations, resulting in a birth rate between 3.3 − 4.1 × 10⁻⁶ yr⁻¹. We discuss the birth rate range obtained in Sect. 5.

The mass of the NS progenitor (M₁) is drawn from Kroupa’s initial mass function (Kroupa 2001), within the range [8–25] M_⊙ (Belczynski & Taam 2008):

$$\begin{array}{c}\hfill P(M_1)\propto M_1^{-2.3}.\end{array}$$(4)

Although our simulation starts with the NS already formed, meaning we are ultimately interested in its initial mass, we still need to know the progenitor mass in order to determine the mass of the companion star. This is because the mass ratio of the binary, $q=\frac{M_2}{M_1}$ (with M₂ being the mass of the secondary star), is drawn from a power-law distribution as suggested by Sana et al. (2012):

$$\begin{array}{c}\hfill P(q)=q^{-0.1}.\end{array}$$(5)

Then, we computed M₂ by deriving q and M₁. Each pulsar has its birth mass, M_NS, i, drawn from a Gaussian distribution such as

$$\begin{array}{c}\hfill P(M_{\mathrm{NS},\mathrm{i}})=\frac{1}{{\sigma }_M\sqrt{2\pi }}{e}^{-{(M_{\mathrm{NS},\mathrm{i}}-{\mu }_M)}^2/(2{\sigma }_M^2)},\end{array}$$(6)

where μ_M = 1.33 M_⊙ and σ_M = 0.09 M_⊙. This distribution is derived from a fit of the Galactic binary NS masses (Özel et al. 2012; Özel & Freire 2016). Our independent sampling of progenitor mass and NS mass removes the statistical correlation whereby more massive progenitors tend, on average, to produce more massive NSs (Sukhbold et al. 2016). This may influence the NS mass distribution obtained in this work. MSPs are very fast rotating NSs. Therefore, for P < 6.5 ms, we can no longer neglect the spin-down caused by gravitational waves (GWs); see Appendix A for more details. We can define a triaxial inertia tensor that is diagonalised on the basis of its principal axes. These axes are denoted by e₁, e₂, e₃ with the corresponding moments of inertia given by I₁ = I, I₂ = I(1 + ϵ₁₂) and I₃ = I(1 + ϵ₁₃) respectively. The terms ϵ₁₂ and ϵ₁₃ characterise the ellipticity, i.e. the degree of deviation from perfect stellar sphericity. However, in our case we assumed we only had biaxial stars, which is sufficient for our purposes, as considering triaxial stars would significantly increase the complexity without providing additional physical insight beyond the key requirement that the stars are deformed. Therefore, we have ϵ₁₂ = 0 and ϵ₁₃ = ϵ ≠ 0. Thus, for each NS, we drew its ellipticity ϵ from a log uniform distribution between [–13,–5] in order to have a value of ϵ between [10⁻¹³, 10⁻⁵]. These limits are suggested by the study of Giliberti & Cambiotti (2022).

The SEVN code also requires the metallicity of the secondary star. As it would necessitate taking into account the environment and the past of the binaries to have an accurate value taken for the metallicity, no assumption was made about it. The metallicity, Z, was taken from a random log-uniform distribution between [0.0002,0.02]. Furthermore, still concerning the secondary star, SEVN also needs the initial ratio Ω_s/Ω_c, which is the spin frequency of the star over its critical spin frequency, we take it from a uniform distribution between [0,1]. We chose a uniform distribution over all possible values of this ratio, as we found that the choice of distribution had little influence on the results.

The binary eccentricity (e) distribution was derived from a power law, as suggested by the study of Sana et al. (2012), where they found this distribution for massive binaries:

$$\begin{array}{c}\hfill P(e)=e^{-0.42}.\end{array}$$(7)

The semi-major axis (a) of the binary was drawn from a distribution that is flat in log(a) for wide binaries and decreases smoothly at close separations, as proposed by Han et al. (2020):

$$\begin{array}{c}\hfill an(a)=\{\begin{array}{cc}{\psi }_{\mathrm{sep}}{\left(\frac{a}{a_0}\right)}^m\hfill & \mathrm{if}a\le a_0\hfill \\ {\psi }_{\mathrm{sep}}\hfill & \mathrm{if}{a}_0<a<a_1,\hfill \end{array}\end{array}$$(8)

where ψ_sep ≈ 0.070, a₀ = 10 R_⊙, a₁ = 5.75 × 10⁶ R_⊙, R_⊙ ≈ 696 340 km and m ≈ 1.2. We do not claim that these distributions represent the exact initial distributions for NS+MS binaries, which are currently unknown. Nevertheless, they cover the relevant parameter space for a and e, and using uniform distributions over the same ranges does not reproduce the observed population as well. Therefore, we retained the distributions of a and e for massive binaries, even though they may not accurately reflect reality. The angle between the rotation axis of the NS and the normal to the orbital plane, which we call α, the spin-orbit angle is assumed to follow an isotropic distribution generated from a uniform distribution U ∈ [0, 1] and given by α = arccos(2 U − 1).

2.4. Neutron star evolution

2.4.1. General evolution model

Because SEVN evolves pulsars in a vacuum as they spin down, we changed the evolution proposed by SEVN with a force-free model evolution taking into account the plasma in the magnetosphere. Further, we adapted the model from Biryukov & Abolmasov (2021) to also take into account the GW emission if the NS reaches a spin frequency Ω ≥ 960 rad.s⁻¹. As soon as the NS reaches Ω < 960 rad.s⁻¹ again, the GWs are neglected. While a fully accurate treatment would require modelling the impact of GWs on the inclination and spin-orbit angle evolution, we neglect these effects as they would greatly complicate the model, while most of the time the GW emission can still be neglected. Thus, in this work the evolution of an NS can be described by the following set of differential equations:

$$\begin{array}{cc}\hfill I\dot{\mathrm{\Omega }}& =n_{1}cos\alpha +n_2+n_3(1+{sin}^2\chi )+n_{\mathrm{GW}}-\dot{I}\mathrm{\Omega },\hfill \end{array}$$(9a)

$$\begin{array}{cc}\hfill I\mathrm{\Omega }\dot{\alpha }& =-n_{1}sin\alpha ,\hfill \end{array}$$(9b)

$$\begin{array}{cc}\hfill I\mathrm{\Omega }\dot{\chi }& =\eta A(\eta ,\alpha ,\chi )n_1{sin}^2\alpha cos\alpha sin\chi cos\chi +n_{3}sin\chi cos\chi ,\hfill \end{array}$$(9c)

where $I=\frac{2}{5}M{R}^2$ is the moment of inertia of the NS (we assume the spherical formula for I, although we consider an ellipticity ϵ for the NS, as this ϵ is very small, between 10⁻¹³ and 10⁻⁵). The torques n₁, n₂, n₃, n_GW are respectively, the averaged torques acting on the NS caused by accretion, magnetic braking due to magnetic field-disc interaction, pulsar’s radiation loss and gravitational braking due to a non-null ellipticity. The coefficient η is set to unity as suggested by Biryukov & Abolmasov (2021). It describes the accretion torque modulation within the spin period. The term A(η, α, χ) is a normalisation function:

$$\begin{array}{c}\hfill A(\eta ,\alpha ,\chi )={[1-\frac{\eta }{2}({sin}^2\chi {sin}^2\alpha +2{cos}^2\chi {cos}^2\alpha )]}^{-1}.\end{array}$$(10)

The torques are defined by

$$\begin{array}{cc}\hfill n_1& ={\dot{M}}_{\mathrm{NS}}{\left(G{M}_{\mathrm{NS}}{r}_{\mathrm{in}}\right)}^{1/2},\mathrm{if}{r}_{\mathrm{in}}<r_{\mathrm{co}},\hfill \end{array}$$(11a)

$$\begin{array}{cc}\hfill n_2& =-\frac{{\mu }^2}{3r_{\mathrm{co}}^3},\mathrm{if}{r}_{\mathrm{in}}<r_{\mathrm{lc}},\hfill \end{array}$$(11b)

$$\begin{array}{cc}\hfill n_3& =-\frac{{\mu }^2}{r_{\mathrm{lc}}^3},\hfill \end{array}$$(11c)

$$\begin{array}{cc}\hfill n_{\mathrm{GW}}& =-\frac{4G{I}^2{\mathrm{\Omega }}^5{ϵ}^2}{5{\pi }^2{c}^5},\mathrm{if}\mathrm{\Omega }\ge 960{\mathrm{rad}.\mathrm{s}}^{-1}.\hfill \end{array}$$(11d)

Here, G is the gravitational constant, M_NS is the mass of the NS, Ṁ_NS is the mass accretion rate of the NS, ${\mu }^2=\frac{4\pi B^2{R}^6}{{\mu }_0}$ is the magnetic moment with the permeability constant μ₀ = 4π × 10⁻⁷ H/m, R = 12 km is the radius of the NS, and c = 3 × 10⁸ m/s is the speed of light in vacuum. The different radii defined, i.e. r_in, r_co and r_lc, are respectively the inner disc radius, the co-rotation radius, and the light cylinder radius:

$$\begin{array}{cc}\hfill r_{\mathrm{in}}& =\xi {\left(\frac{{\mu }^4}{2G{M}_{\mathrm{NS}}{\dot{M}}_{\mathrm{NS}}^2}\right)}^{1/7},\hfill \end{array}$$(12a)

$$\begin{array}{cc}\hfill r_{\mathrm{co}}& ={\left(\frac{G{M}_{\mathrm{NS}}}{{\mathrm{\Omega }}^2}\right)}^{1/3},\hfill \end{array}$$(12b)

$$\begin{array}{cc}\hfill r_{\mathrm{lc}}& =\frac{c}{\mathrm{\Omega }},\hfill \end{array}$$(12c)

where ξ ∼ 0.5 is a dimensionless factor accounting for the detailed disc structure (Bessolaz et al. 2008).

It is not straightforward to visualise the temporal evolution of a pulsar in the P−Ṗ diagram from this system of differential equations. Therefore, we provide in Sect. 2.7 an example of the typical evolutionary track of a pulsar that becomes a recycled pulsar.

2.4.2. Spin-down

When the NS is not accreting matter, n₁ and n₂ are set to zero, and the star spins down as a result. Additionally, n_GW is also set to zero when Ω < 960 rad.s⁻¹. When these conditions are satisfied, we can apply the simplified set of equations from Sautron et al. (2024), which consider only the pulsar’s radiative losses. This approach allows for faster computation compared to solving the full system of Eqs. (9a), (9b), and (9c). The equations are given by

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

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

where the quantities with a subscript of 0 indicate their value at birth. Those without a subscript display their current value at present time. The time t represents the true age of the pulsar and K = 4πR⁶/Iμ₀c³. Equation (13a) is the integral of motion between Ω and χ, the inclination angle, and Eq. (13b) allowed us to find the evolution of χ in a dipolar configuration with a magnetic field decay. This system was found by Philippov et al. (2014) for a spherically symmetric NS with a constant magnetic field and Dirson et al. (2022) adapted Eq. (13b) for a decaying prescription. Finally, when there is no accretion and Ω ≥ 960 rad.s⁻¹, we solve still Eqs. (9a) and (9c), with n₁ = n₂ = 0.

2.4.3. Spin-up

As mentioned earlier, the NS can spin-up only if it is accreting matter, as n₁ is the only torque that is positive. Therefore, only n₁ can increase Ω. Two conditions are necessary for the NS to accrete matter. First, the companion must reach its Roche lobe, and then the NS needs to have its Keplerian angular velocity at the inner disc radius (Ω_K|rin) be greater than its co-rotation angular velocity (Ω_co):

$$\begin{array}{c}\hfill {\mathrm{\Omega }}_{\mathrm{diff}}={\mathrm{\Omega }}_{\mathrm{K}|r_{\mathrm{in}}}-{\mathrm{\Omega }}_{\mathrm{co}}.\end{array}$$(14)

Thus, Ω_diff > 0 is necessary for the NS to be accreting matter, the inner disc radius r_in is defined by Eq. (12a). When matter reaches the inner disc radius, the magnetic pressure dominates, forcing the material to follow the magnetic field lines and funnel onto the NS’s polar caps. However, if Ω_diff < 0, the NS enters the propeller regime: the magnetic field lines at the magnetic radius rotate faster than the local Keplerian velocity, preventing accretion and expelling the matter away from the star (Abolmasov et al. 2024).

2.4.4. Magnetic field evolution

The evolution of NS magnetic fields remains a subject of debate. Some early pulsar population synthesis studies suggested that the magnetic field does not evolve over time (Bhattacharya et al. 1992; Faucher-Giguère & Kaspi 2006). However, several physical mechanisms support the idea of magnetic field decay, including ohmic dissipation in the crust, ambipolar diffusion in the core, crustal vortex motion, internal turbulence, and energy loss through radiation (Srinivasan et al. 1990; Cumming et al. 2004; Dall’Osso et al. 2012; Viganò et al. 2013). A recent study that includes ohmic heating (Igoshev & Popov 2024), along with period derivative measurements and spectral data from the NS RX J0720.4-3125, further suggests a rapid decay of its magnetic field. Moreover, more recent population synthesis studies have shown that models incorporating magnetic field decay provide a great match to the observed pulsar population (Gullón et al. 2014; Sautron et al. 2024). Therefore, we adopted a magnetic field decay following a power-law prescription:

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

where α_d is a constant parameter controlling the speed of the magnetic field decay, B₀ the initial magnetic field, B_min the minimum magnetic field reachable by the NS, t corresponds to the elapsed time since the birth of the NS and τ_d the typical decay timescale. We use α_d = 1.5 (as in Sautron et al. 2024 for the normal pulsars). Concerning τ_d, we considered three different values randomly chosen for the generated pulsar, 0.05, 1.4, and 2.3 Myr, with a probability of 0.77, 0.115, and 0.115 to have these values, respectively. These three probabilities effectively approximate the probability density function of τ_d and were sufficient to reproduce the recycled pulsar population. Determining the full distribution of τ_d would require further investigation beyond the scope of this study. Physically, the τ_d distributions represent the range of possible evolutionary trajectories in the P–Ṗ diagram for recycled pulsars, similar to Fig. 10 in Viganò et al. (2013). Thus, each pulsar in the simulation follows one of three possible evolutionary paths, an empirical approach that reproduces the P–Ṗ distribution better than using a single path. These discrete values were selected purely by trial and error. Whether the NS is accreting or not, the magnetic field is decaying in this model; however accretion-induced field decay was also taken into account.

As the NS accretes matter from its companion, this matter goes along the magnetic field lines to the polar caps. Progressively, the matter accreted can apply enough pressure to compress and bury the magnetic field under the NS surface (Payne & Melatos 2007; Wang et al. 2012). This accretion-induced field decay process can be described by

$$\begin{array}{c}\hfill B_{\mathrm{AI}}(t)=(B(t)-B_{\min })e^{-\frac{\mathrm{\Delta }{M}_{\mathrm{NS}}}{\mathrm{\Delta }{M}_d}}+B_{\min },\end{array}$$(16)

where B(t) is computed via Eq. (15) to model the decay caused by the other physical mechanisms mentioned earlier, ΔM_NS is the amount of accreted mass by the NS and ΔM_d is the magnetic field decay mass scale. Typically the accreted mass by an NS is less than 0.5 M_⊙ (Özel & Freire 2016; You et al. 2025) and we adopt ΔM_d = 0.025 M_⊙ as in the work of Chattopadhyay et al. (2020) which is the optimal value found in their study.

To prevent the magnetic field of NSs from dropping to unrealistically low values, we imposed for each pulsar a minimum magnetic field value, B_min, randomly drawn from a uniform distribution in the range [10³, 10⁴] T. This range is motivated by the study of Zhang & Kojima (2006), who investigated the lower limits that NS magnetic fields can realistically attain.

2.5. Spatial evolution

The motion of each binary’s centre of mass within the Galactic potential is evolved in the gravitational potential Φ. They are subject to an acceleration $\ddot{\mathit{x}}$ given by

$$\begin{array}{c}\hfill \ddot{\mathit{x}}=-\mathbf{\nabla }\mathbf{\Phi }.\end{array}$$(17)

We integrate numerically Eq. (17) with a position extended Forest Ruth-like (PEFRL) algorithm (Omelyan et al. 2002), a fourth order integration scheme. Concerning the galactic potential, the Galaxy is divided in four distinct regions, which results in four different potentials: the bulge Φ_b and the disc Φ_d which have the form proposed by Miyamoto & Nagai (1975), the dark matter halo Φ_h which has the form proposed by Navarro et al. (1997) and the nucleus Φ_n represented by a Keplerian potential. The total potential of the Milky Way Φ_tot is the sum of these potentials:

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathit{tot}}={\mathrm{\Phi }}_b+{\mathrm{\Phi }}_d+{\mathrm{\Phi }}_h+{\mathrm{\Phi }}_n.\end{array}$$(18)

See Sautron et al. (2024) for a detailed analysis of the PEFRL scheme, its accuracy, and for the detailed expression of the galactic potentials used.

Although SEVN does not integrate the full orbital motion of binaries with respect to their centre of mass, it computes the evolution of the semi-major axis and eccentricity using prescriptions for mass exchange within the system, following Eqs. (15) and (16) of Hurley et al. (2002). The dominant processes affecting these orbital parameters are tides, which are modelled through Eqs. (34)–(36) of Iorio et al. (2023, see Sects. 2.3.1 and 2.3.4 of their paper for more details about wind mass transfer and tides respectively). In addition, SEVN accounts for the evolution of the semi-major axis and eccentricity driven by GW emission. Peters (1964) describes in his work the orbital decay and circularisation caused by GWs with the following equations:

$$\begin{array}{cc}\hfill \frac{\mathit{da}}{\mathit{dt}}& =-\frac{64G^3{m}_1{m}_2(m_1+m_2)}{5c^5{a}^3{(1-e^2)}^{7/2}}(1+\frac{73}{24}{e}^2+\frac{37}{96}{e}^4),\hfill \end{array}$$(19a)

$$\begin{array}{cc}\hfill \frac{\mathit{de}}{\mathit{dt}}& =-\frac{304G^3{m}_1{m}_2(m_1+m_2)}{15c^5{a}^3{(1-e^2)}^{5/2}}(1+\frac{121}{304}{e}^2),\hfill \end{array}$$(19b)

where G is the gravitational constant, m₁ and m₂ are the masses of the two stars in the binary, c is the speed of light, a is the semi-major axis and e is the eccentricity of the orbit. In SEVN, these equations are solved when the GW merger timescale, t_merge becomes shorter than the Hubble time. See Appendix D of Iorio et al. (2023) for more details about the computation of the GW merger timescale.

2.6. Death valley

Spread in the parameter values of the model causes significant variations in the death line, and thus a death valley rather than a single death line describes the condition for pulsar extinction in the P−Ṗ diagram. Nonetheless, this death valley must be based on a death line definition, we choose to extend the death line from Mitra et al. (2020), originally defined for normal pulsars to access the criteria for generation of coherent radio emission given as follows:

$$\begin{array}{c}\hfill {\dot{P}}_{\mathrm{line}}=\frac{3.16\times {10}^{-19}{T}_6^4{P}^2}{{\eta }^{2}b{cos}^2{\chi }_l}·\end{array}$$(20)

Here, T₆ = T/10⁶ K, where T corresponds to the surface temperature of the polar cap and P is the spin period of the pulsar in s. The parameter η = 1 − ρ_i/ρ_GJ is the electric potential screening factor due to the ion flow, where ρ_i corresponds to the ion charge density and ρ_GJ to the Goldreich-Julian charge density above the polar cap. The quantity b is the ratio of the actual surface magnetic field to the dipolar surface magnetic field and χ_l is the angle between the local magnetic field and the rotation axis. In order to obtain the death valley in Fig. 1, the parameters described above are drawn from distributions, found and described in Sautron et al. (2024).

Fig. 1. Diagram of P−Ṗ of the observed pulsars considered in this work along with the death line (solid green line) and the death valley (shaded green area).

2.7. Evolution in the P−Ṗ diagram of a recycled pulsar

To better visualise how a pulsar evolves in the P−Ṗ diagram to become a recycled pulsar, Fig. 2 shows the track followed by a pulsar evolving in the simulation, starting with an initial period P₀ = 9.7 × 10⁻² s, period derivative Ṗ₀ = 2.7 × 10⁻¹³ s/s, and mass M_NS, 0 = 1.38 M_⊙, its companion being a ∼8 M_⊙ main sequence star. Since the mass ratio of these two stars, $q=\frac{M_2}{M_{\mathrm{NS}}}=5.8$, exceeds the critical mass ratio for main sequence stars in the simulation, q_c = 3.0, it means that the mass transfer would be unstable. This implies that a CE phase must have occurred, allowing the secondary star to lose enough mass in the envelope, which was then ejected, for the system to later enter a stable mass-transfer phase. For further details on binary evolution, see Tauris & van den Heuvel (2023). We can basically distinguish three major phases in the evolution of this pulsar: first, the pulsar was dominated by the spin-down (caused mostly by the torque n₃ presented in 2.4) and evolved towards a spin period close to 10 s. Then, when the accretion started (after 0.21 Gyr of evolution) we can clearly see that the pulsar, which passed the death line before, comes back from the commonly called graveyard of pulsars to have a very rapid rotation, close to the millisecond, thanks to the spin-up (caused by the torque n₁ presented in 2.4) that dominates this phase. Finally, when the accretion ends (the accretion lasted ∼0.33 Gyr), the pulsar can continue to evolve by being dominated once again by spin-down. However, compared to when the pulsar was young, the spin-down is much weaker, due to the reduced value of the magnetic field, which is now much lower than at its birth, that is why ∼10 Gyr after the end of the accretion phase, the pulsar barely moved in the P−Ṗ diagram. The pulsar finished its evolution with P = 1.9 × 10⁻³ s, Ṗ = 1.1 × 10⁻²⁰ s/s, M_NS = 2.3 M_⊙ (thus, the pulsar accreted ∼0.9 M_⊙) and ended up isolated as its companion was ejected after its supernova explosion.

Fig. 2. Evolution in the P−Ṗ diagram of a pulsar in the simulation that becomes an MSP.

3. Detection model

In this section we address the detectability of our simulated pulsars in both radio and γ-ray wavelengths. First of all, in this work we exclusively consider two radio surveys, the FAST Galactic Plane Pulsar Snapshot (GPPS) survey (Han et al. 2021) and the Parkes Multibeam Pulsar Survey (PMPS; Manchester et al. 2001), as these surveys detected the most MSPs. For the γ-ray wavelengths, we are interested in the comparison of the results of our simulation with the detections obtained by the Fermi/LAT instrument. Therefore, we compared the γ-ray pulsars obtained in our simulation with the pulsars from the Third Fermi Large Area Telescope Catalog (3PC).

For each pulsar, we verified whether it met the detection criteria based on three factors. The first is the beaming fraction, which represents the portion of the sky swept by the radiation beam and depends on the wavelength (radio or γ-ray), the spin rate, geometry, and emission region. Before detailing how the radio and γ-ray beaming fractions are computed, we first define the relevant angles.

The angle between the line of sight and the rotation axis is $\zeta =(\widehat{{\mathit{n}}_{\mathbf{\Omega }},\mathit{n}})$, where n_Ω = Ω/||Ω|| is the unit vector along the rotation axis, and n is the unit vector along the line of sight. The inclination angle between the rotation axis and the magnetic moment is $\chi =(\widehat{{\mathit{n}}_{\mathbf{\Omega }},\mathit{\mu }})$, with μ the unit vector along the magnetic moment. Finally, the impact angle $\beta =(\widehat{\mathit{\mu },\mathit{n}})$ is related by χ + β = ζ.

We chose an isotropic distribution for the Earth viewing angle, ζ, as well as for the orientation of the unitary rotation vector. The Cartesian coordinates of the unit rotation vector, n_Ω, are (sin θ_nΩ cos ϕ_nΩ, sin θ_nΩ sin ϕ_nΩ, cos θ_nΩ). We set the Sun’s position at (x_⊙, y_⊙, z_⊙) = (0 kpc, 8.5 kpc, 15 pc) (Siegert 2019). The coordinates for n are

$$\begin{array}{c}\hfill \mathit{n}=(\frac{x-x_{\odot }}{d},\frac{y-y_{\odot }}{d},\frac{z-z_{\odot }}{d}).\end{array}$$(21)

To compute the pulsar distance from Earth, we used the formula for the distance,

$$\begin{array}{c}\hfill d=\sqrt{{(x-x_{\odot })}^2+{(y-y_{\odot })}^2+{(z-z_{\odot })}^2}.\end{array}$$(22)

3.1. Radio detection model

The beaming fraction in radio depends on the half opening angle of the radio emission cone, ρ. Since our selected sample of observed pulsars includes a few pulsars with P > 0.1 s, the half opening angle, ρ, of simulated pulsars that reach a similar value is computed according to Lorimer & Kramer (2004) by

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

where h_em is the emission height, P is the pulsar spin period, and c is the speed of light. The emission height was assumed to be constant with an average value of h_em = 3 × 10⁵ m, based on observational estimates from a sample of pulsars (Weltevrede & Johnston 2008; Mitra 2017; Johnston & Karastergiou 2019; Johnston et al. 2023). The cone half-opening angle, ρ, as defined in Eq. (23), applies only to the last open field lines of a dipolar magnetic field and is valid primarily for slow pulsars, where the emission altitude is high enough for higher-order multi-polar components to be negligible. For fast pulsars, however, relativistic effects due to rapid rotation become significant and must be taken into account (Xilouris et al. 1998). Furthermore, little is known about the characteristics of the half-opening angle and the emission height in fast pulsars. For this reason, and since we adopted a purely dipolar geometry for all pulsars, even though multi-polar components are expected to play a role in reality, we retained a formula similar to Eq. (23) following the prescription of Kijak & Gil (1998, 2003) for ρ:

$$\begin{array}{cc}\hfill \rho & =1.{24}^{°}{r}_{\mathrm{alt}}^{1/2}{P}^{-1/2},\hfill \end{array}$$(24a)

$$\begin{array}{cc}\hfill f(r_{\mathrm{alt}})& =\frac{1}{{\sigma }_{\mathrm{alt}}\sqrt{2\pi }}{e}^{-{(r_{\mathrm{alt}}-{\mu }_{\mathrm{alt}})}^2/\left(2{\sigma }_{\mathrm{alt}}^2\right)},\hfill \end{array}$$(24b)

where r_alt denotes the emission altitude in units of NS radii, and μ_alt and σ_alt are respectively the mean and standard deviation of the Gaussian distribution (f(r_alt)) from which r_alt is drawn. Running multiple simulations showed that adopting μ_alt = 8 R_NS and σ_alt = 4 R_NS provides the best agreement with observations, meaning that both radio pulse profiles and the fraction of γ-ray pulsars that could be visible in radio are consistent with the observations.

A pulsar can be detected in radio if β = |ζ − χ|≤ρ or β = |ζ − (π − χ)| ≤ ρ, corresponding to visibility from the northern or southern magnetic hemisphere, respectively. An additional useful quantity is the radio pulse width, w_r, which is computed following Lorimer & Kramer (2004) as

$$\begin{array}{c}\hfill cos(\rho )=cos(\chi )cos(\zeta )+sin(\chi )sin(\zeta )cos(w_r/2).\end{array}$$(25)

For radio flux density, we use a formula similar to that in Johnston et al. (2020), designed for 1.4 GHz observations by the Parkes telescope in the southern Galactic plane (Kramer et al. 2003; Lorimer et al. 2006; Cameron et al. 2020) and Arecibo in the northern plane (Cordes et al. 2006). The main difference in this work lies in the power-law dependence on the spin-down luminosity $\dot{E}=-I\mathrm{\Omega }\dot{\mathrm{\Omega }}$. Johnston & Karastergiou (2017) found that using Ė^0.5 overpredicted young pulsars, while Ė⁰ overrepresented old ones; they adopted Ė^0.25 as a compromise. In our case, we found that lowering the exponent to 0.125 improved agreement with observations. Thus, the flux density formula is

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

where d is the distance in kpc and F_j is the scatter term which is modelled as a Gaussian with a mean of 0.0 and a variance of σ = 0.2. The detection threshold in radio is set by the signal-to-noise ratio defined by

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

We directly compute the radio flux, without computing the luminosity in radio that overlooks the fact that the luminosity received will depend on the geometry of the beam. A pulsar is detected in radio if its signal-to-noise ratio exceeds the threshold imposed by FAST GPPS Han et al. (2021) or PMPS Manchester et al. (2001). Indeed, in this work we focus exclusively on these two radio surveys. The region of observations of these surveys are taken into account; see Table 1. S/N depends on $S_{\mathrm{survey}}^{\min }$ which is the minimum detectable flux, which is related to the instrumental sensitivity, the period of rotation P and the observed pulse profile width of radio emission ${\stackrel{\sim }{w}}_r$. To compute the observed pulse profile width, we adopted the same formula as in Cordes & McLaughlin (2003):

$$\begin{array}{cc}\hfill w_r& =\sqrt{{(w_{r}P/2\pi )}^2+{\tau }_{\mathrm{samp}}^2+{\tau }_{\mathrm{DM}}^2+{\tau }_{\mathrm{scat}}^2},\hfill \end{array}$$(28a)

$$\begin{array}{cc}\hfill S_{\mathrm{survey}}^{\min }& =\frac{C_{\mathrm{thres}}[T_{\mathrm{sys}}+T_{\mathrm{sky}}(l,b)]}{G\sqrt{N_{p}Bt}}\sqrt{\frac{{\stackrel{\sim }{w}}_r}{P-{\stackrel{\sim }{w}}_r}}·\hfill \end{array}$$(28b)

Equation (28a) accounts for interstellar medium (ISM) density fluctuations, including dispersion (τ_DM) and scattering (τ_scat) caused by interactions with free electrons in the Galaxy during the radio pulse’s propagation. Instrumental effects are also considered through the sampling time, τ_samp. More details on these parameters are provided in Appendix B.

In Eq. (28b), where we detail the expression of the minimum detectable flux, C_thres is the detection threshold of the survey, T_sys is the system temperature in K, G is the telescope gain in K.Jy⁻¹, N_p is the number of polarisations, B is the total bandwidth in Hz and t is the integration time in s. The quantities which have just been presented, are the survey characteristics, which depend on whether we consider FAST GPPS or PMPS (see Table 1). T_sky is the sky background temperature in K, at the longitude l and latitude b. Remazeilles et al. (2015) provided a refined version of the temperature map of Haslam et al. (1982). Since the data were obtained at 408 MHz, in order to rescale to the correct observing frequencies, a power-law dependence is assumed with the form T_sky ∝ f^−2.6 (Johnston et al. 1992). We use Eq. (28b) to compute the minimum flux required for radio detection.

3.2. γ-ray detection model

The γ-ray emission model is based on the striped wind scenario, where γ-ray photons are produced in the current sheet of the pulsar wind. Detection requires the observer’s line of sight to lie near the equatorial plane, with the inclination angle satisfying |ζ − π/2|≤χ. The γ-ray luminosity is taken from the study by Kalapotharakos et al. (2019), who found it follows a fundamental plane relation. Their full 3D model depends on the magnetic field B, the spin-down luminosity Ė, and the cut-off energy ϵ_cut. Since our PPS does not include ϵ_cut, we adopted the simplified 2D version of their model:

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

Here, Ė is the spin-down luminosity, and the associated γ-ray flux detected on Earth is computed with

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

where f_Ω is a factor depending on the emission model reflecting the anisotropy. For the striped wind model, Pétri (2011) showed that this factor varies between 0.22 and 1.90. In this work, we use an updated version of Fig. 7 of Pétri (2011), and it allowed us to obtain f_Ω when the viewing angle, ζ, and the inclination angle, χ, are known. Therefore, this update allowed us to get more realistic values for f_Ω between 0.13 and 316. See Appendix C for more details on using the force-free model. The pulsar is detected in gamma depending on the instrumental sensitivity as described below.

Fig. 3. Diagram of P−Ṗ for observed and simulated pulsars in different wavelengths. Left panel: P−Ṗ diagram for the detected pulsars in radio surveys only in the simulations in red along with the observations in radio for the FAST GPPS or PMPS survey in blue. Middle panel: P−Ṗ diagram for the detected pulsars in γ-ray surveys only in the simulations in green, along with the observed pulsars in γ-ray from Fermi surveys only in blue. Right panel: P−Ṗ diagram for the simulated pulsars detected both in radio and γ-ray surveys in purple, along with observations for the FAST GPPS or PMPS survey simultaneously with the γ-ray observations from Fermi in blue. For all the panels, the death line is represented as a green solid line, and the death valley is the shaded green area.

The γ-ray sensitivity is determined based on the expected performance of the Fermi/LAT instrument³. For each simulated pulsar, we consult the all-sky sensitivity map provided by Fermi/LAT⁴. At the position defined by the Galactic longitude (l) and latitude (b) of the pulsar, we retrieve F_min, the minimum detectable flux at that location. If Fermi/LAT did not observe the region, this is also taken into account. A pulsar is considered detectable in γ-rays if its γ-ray flux (F_γ) exceeds the local threshold (F_min).

4. Results

4.1. Radio versus γ-ray pulsars

The first goal of this study was to match the P−Ṗ diagram between simulation and observations in order to assess whether the recycling process could primarily explain the observed pulsar population of mildly and fully recycled pulsars in the P−Ṗ diagram (i.e. pulsars with B ≤ 5 × 10⁶ T). By using the parameters for the various distributions at birth and the evolution model presented in Sect. 2, we obtained the three P−Ṗ diagram in Fig. 3. The left panel displays the diagram for pulsars detected only in radio surveys, the middle panel for those detected only in γ-rays surveys, and the right panel for pulsars detected in both radio and γ-rays surveys. In each case, observations are shown alongside simulations. For the radio-only population (left panel), the Kolmogorov–Smirnov (KS) test (Smirnov 1948) yields a p-value of 5.6 × 10⁻⁵ for the spin period P distribution and 0.2 for the spin-down rate Ṗ distribution. This suggests that the null hypothesis is rejected for P but not for Ṗ. Nevertheless, the simulated population matches observations well for P < 0.1 s, with p-values greater than 0.05 for both P and Ṗ, indicating that few recycled pulsars are observed above this value, meaning this channel of formation only allows one to obtain a few mildly recycled pulsars, but not all of them. In the ATNF catalog, 167 radio pulsars detected by the FAST GPPS or PMPS surveys have measured values of both P and Ṗ. Additional detections exist but are sometimes missing a Ṗ measurement – such detections were excluded from the comparison with our simulation. Other additional detections may also exist that are not yet included in the catalog. A more complete list of sources can be found online⁵. In our simulation, we detect 176 radio pulsars, which agrees reasonably well with the number of radio pulsars discovered by FAST GPPS and PMPS. Furthermore, the death valley eliminates ∼100 pulsars, and most of them were well below the valley. Regarding, the middle panel (pulsars in γ-ray surveys only), most of the observed pulsars appear to be reproduced in the simulation, though there is a noticeable over-detection for Ṗ > 3 × 10⁻²⁰. From 3PC, 137 pulsars were selected for comparison (still with the criterion: B ≤ 5 × 10⁶ T), of which 20 were also detected by the selected radio surveys, leaving 117 as pulsars in γ-ray surveys only. In the simulation, 98 such pulsars were detected with Ṗ < 3 × 10⁻²⁰. When comparing this observed sample to the simulated one (restricted to Ṗ < 3 × 10⁻²⁰), the KS test yields p-values of 0.15 for Ṗ and 0.25 for P, indicating good agreement. Furthermore, the number of detected pulsars in both datasets is very similar. The overabundance of simulated pulsars with Ṗ > 3 × 10⁻²⁰ could correspond to unidentified sources in the Fourth Fermi-LAT catalogue of γ-ray sources (4FGL), or may reflect limitations of the striped wind model. These interpretations are further discussed in Sect. 5. Finally, the population of pulsars detected in both radio and γ-ray surveys (right panel) is too small for a robust KS test. However, the simulated pulsars occupy a similar region of the P–Ṗ space as the observed ones. The simulation does show an excess of approximately 20 pulsars compared to observations, which corresponds well to the number of radio MSPs associated with γ-ray sources that have not yet been observed to pulsate in the 4FGL catalogue (Abdollahi et al. 2022). These observed sources may simply not have shown detectable pulsations so far, and further investigations could reveal their MSP nature, which would be consistent with our simulation results.

Although the middle panel displays what we classified as pulsars in γ-ray surveys only, we still computed their expected radio flux using Eq. (26). This choice is motivated by the fact that, once a pulsar is detected in γ rays, it is often followed up with dedicated radio observations. As a result, most of the γ-ray pulsars eventually have some form of radio measurement associated with them. However, these follow-up observations are not linked to any specific survey. As a result, in the ATNF data, many γ-ray pulsars that do emit in radio are not listed as being associated with a radio survey, which is the criterion we used to classify pulsars as either radio-loud or radio-quiet in Fig. 3. Consequently, ∼71% of the objects we label as pulsars in γ-ray surveys only in the simulation actually have radio fluxes greater than 30 μJy and a geometry that would allow one to see pulsations, i.e |ζ − χ|≤ρ, which in real observations would almost always lead to a radio detection (Smith et al. 2023). Indeed, in the sample of 117 pulsars used in this work from the 3PC catalogue labelled as pulsars in γ-ray surveys only, only 6 of them actually lack any radio detection, either from surveys or from targeted observations. This corresponds to 94.9% of γ-ray pulsars not seen in a radio survey, which have a radio counterpart.

We compare the simulated and observed γ-ray light-curve peak separation, Δ, in Fig. 4. The separation, Δ, was computed according to Pétri (2011) with

$$\begin{array}{c}\hfill cos\left(\pi \mathrm{\Delta }\right)=|cot\zeta cot\chi |.\end{array}$$(31)

Fig. 4. Distribution of the γ-ray peak separation of the simulated population along with the observations from the 3PC catalogue.

Here, ζ is the angle between the line of sight and the rotation axis, and χ is the magnetic inclination angle. All distributions are normalised to yield probability distribution functions (PDFs) by dividing each PDF by the total number of pulsars–observed or simulated. Since the signal is periodic and the definition of the main peak is arbitrary, the peak separation (Δ) is restricted to the interval [0, 0.5]; values larger than 0.5 are mapped to 1 − Δ via a phase shift. Figure 4 shows a striking agreement between the observed and simulated γ-ray peak separation distributions, indicating that the striped wind model reproduces this geometric feature well, as was also the case for the normal pulsars in Sautron et al. (2024). Moreover, for the observed sample, separations smaller than 0.15 are likely unresolved, which explains why this region shows poorer agreement with the observations. Nevertheless, the overall trend of the curve remains consistent between the observations and the simulations.

The similarity in emission geometry between the simulation and the observations is also evident in Fig. 5, where we compare the width of the radio profile as a function of the spin period, P, for the simulated pulsars above detection threshold and the observed ones. A linear fit in log scale to both datasets yielded consistent trends, with w_r ∝ P^{−0.34 ± 0.04} for the simulation and $w_{10}^{\mathrm{ATNF}}\propto P^{-0.33\pm 0.05}$ for the observations. The simulation successfully reproduces the full range of observed pulse widths, further indicating that the adopted radio emission geometry is realistic. Although one might expect a larger ratio between the two curves (here ∼0.97) if the radio profiles were purely Gaussian, since we compare the full simulated profile width according to Eq. (25) with the observed width measured at 10% of the peak intensity, a ratio closer to ∼1.7 would be anticipated. Nevertheless, this difference can likely be attributed to the use of a dipolar geometry, to uncertainties in the measurement of $w_{10}^{\mathrm{ATNF}}$, and to the skewness of some observed profiles.

Fig. 5. Width of the radio profile plotted against the spin period, P, of the radio-detected pulsars in the simulation along with the observations.

4.2. Companion

Since every pulsar in the simulation begins with a companion, it is instructive to investigate their eventual fate. The left panel of Fig. 6 shows the types of companions associated with NSs at the end of the simulation (noting that some NSs end up isolated). Most of the detected simulated systems end as isolated NSs or with a helium WD (HeWD) or carbon–oxygen WD (COWD) companion. Only a few systems retain other types of companions, such as non-degenerate stars, other NSs, oxygen–neon WDs (ONeWDs) or even a BH in some cases (although this simulation did not produce any BH companion). The right panel displays the initial companion mass as a function of the accretion duration. This plot informs us quite clearly that above an initial mass of 8 M_⊙ for the companion, the survival of the binary is very rare. Most of the pulsars that end up isolated, had a companion whose initial mass was above 8 M_⊙. Therefore, because of the supernova explosion, most of the time the binary is disrupted by the kick velocity of the companion at the remnant birth. This explains the rarity of NS-NS or NS-BH systems: their formation requires a companion with an initial mass above 8 M_⊙, followed by a sufficiently low natal kick at the NS’s birth to prevent the binary from being disrupted.

Fig. 6. Companion outcomes in the pulsar population. Left panel: Distribution of the type of the pulsars companions. HeWD stands for helium white dwarf, ‘isolated’ means the pulsar does not have a companion anymore, COWD stands for carbon-oxygen white dwarf, NotARemnant means the companion is a non-degenerate star, NS CCSN indicates an NS born after a core-collapse supernova, NS ECSN refers to an NS born after an electron capture supernova, and ONeWD stands for oxygen-neon white dwarf. Right panel: Initial mass of the pulsars’ companions plotted against the duration of the accretion phase, with the final type of the companion displayed. MS stands for main sequence star, and NS stands for NS.

A further noticeable aspect of Fig. 6 is the fact that the initial mass of the companion is not correlated with the duration of the accretion phase of the binary. Basically, most of the binaries have an accretion phase which lasts between 10⁸ and 10⁹ yr. Furthermore, with Eq. (9b) we deduced a typical timescale of alignment for the rotation axis of the pulsar and the normal to the orbital plane, $t_{\mathrm{align}}^{\mathrm{orbit}}$:

$$\begin{array}{c}\hfill t_{\mathrm{align}}^{\mathrm{orbit}}=\frac{I\mathrm{\Omega }}{n_1},\end{array}$$(32)

where n₁ is the accretion torque, I is the moment of inertia, and Ω is the frequency of rotation of the pulsar. If we compute this timescale with M_NS = 1.4 M_⊙, Ṁ_NS = 10⁻⁹ M_⊙.yr⁻¹, R = 12 km, B = 5 × 10⁶ T, Ω = 1 rad.s⁻¹, and I = 10³⁸ kg.m², which are typical values for a pulsar that has been spinning down for a long time without accretion, it gives $t_{\mathrm{align}}^{\mathrm{orbit}}\sim 5.9$ kyr. Even if we change these values, in many cases, the typical duration of spin-orbit alignment is going to be lower than the typical duration of accretion, meaning most of the population has its spin-orbit aligned. Indeed, we already showed these results in Lorange et al. (2026, see Fig. 7), corresponding to ∼80% of the population with the spin-orbit angle, α, that is less than 10°. This supported the possibility of detecting pulsars with a misaligned spin-orbit (even though simulations show that such cases are relatively rare), as this study specifically investigated the spin-orbit alignment hypothesis.

4.3. Distribution of NS masses

Studying the mass that NSs can reach is very important in order to constrain the equation of state (EoS; Özel & Freire 2016). In this work we are able to follow the mass evolution of the NSs, and especially look at what is the distribution of mass characterising the recycled pulsars detected in the simulation. In Fig. 7, the PDF of the detected pulsar mass obtained in the simulation is represented. Linares (2020) gives a review of the NS mass distribution observed. An histogram including 44 mass measurements of pulsars was shown (Fig. 4 of Linares 2020), and at that time it was very new to obtain measurements of pulsars with masses greater than 2 M_⊙, and for most of them they were spiders system. To this day, the number of measurements of mass of pulsars did not increase much⁶. Nevertheless, many of the observed mass measurements correspond to young pulsars, which clearly did not have time to undergo recycling. Moreover, the total number of available measurements remains limited. In contrast, in this type of simulation, we can obtain the mass of each recycled pulsar. The results suggest that due to accretion, the mass distribution of recycled pulsar is centred around 1.8 M_⊙, with pulsars potentially reaching up to 2.7 M_⊙. This is consistent with recent mass measurements of spider pulsars mentioned earlier. Therefore, the simulation indicates that if more masses of recycled NSs were measured, we could expect to find significantly more high masses compared to normal pulsars. Although the maximum NS mass imposed in SEVN was 3 M_⊙, no pulsar in the simulation reached this value. There is also the possibility that the pulsars in the simulation accrete too much matter, as we did not observe many high-mass pulsars. Therefore, to diminish the accretion in the simulation, it would require larger τ_d so that the magnetic field could decay more slowly. This, in turn, would increase the inner disc radius r_in for a longer period, thereby shortening the accretion phase. We tested this hypothesis in the simulation, but it failed to reproduce a realistic P−Ṗ diagram, making this scenario unlikely. Besides, as already mentioned, the number of measured masses of NSs is small (< 100). Overall, this study supports the recent high mass measurements of spider pulsars and suggests that the maximum mass of NSs may be higher than the 2.28 M_⊙ proposed by Ruiz et al. (2018); otherwise, the model would not allow the P−Ṗ diagram to be reproduced so well.

Fig. 7. Histogram of the final mass of the recycled pulsars in the simulation. The vertical dashed black line represents the theoretical maximum mass reachable by an NS from the study of Ruiz et al. (2018).

4.4. Spatial distribution

Very interesting result concerns Fig. 8, where we can see the distributions in Galactic coordinates for the detected pulsars in the simulations and in the observations. It is striking that we recover well the positions of the observed pulsars with our simulations. However, we detected more pulsars especially at latitudes between −15° and 15° in the simulation. Therefore, we could expect many unidentified sources of 4FGL to lie at these latitudes. Furthermore, regarding the measured positions, there are uncertainties which can be greater than 20% on their measurement (Verbiest et al. 2012; Igoshev et al. 2016; Yao et al. 2017), which one must have in mind while comparing the simulated and observed positions of pulsars. In addition, this simulation can give us an estimation of the number of recycled pulsars born in the Galactic field which end up in the Galactic centre. Two majors hypothesis exist regarding the origin of the GeV excess. Either it is of a self-annihilating dark matter particle origin (Calore et al. 2015; Ackermann et al. 2017; Hooper 2022), or it is of pulsar origin (Yang & Aharonian 2016; Malyshev 2024; Sautron et al. 2024). Moreover Yang & Aharonian (2016) computed a total luminosity of 10³⁰ W on the inner 10° ×10° region of the Galactic centre, which corresponds to a projected square of 1.4 kpc × 1.4 kpc around the Galactic centre (assuming a distance of 8.5 kpc from the Sun), i.e. a radial extent of ∼0.7 kpc. In the simulations, it would correspond to a total of ∼4.6 × 10²⁸ W generated by ∼190 γ-ray pulsars on average in a 0.7 kpc radius from the Galactic centre (whether these pulsars are detected or not, and usually there is between 0 and 2 γ-ray detection in that area in the simulations). Therefore, the simulation hint to the fact that MSPs that were originally born in the Galactic field seem to only contribute to ∼5% of the GeV excess. In addition, Sautron et al. (2024) also found that young γ-ray pulsars could contribute to this excess with their detection prospects, and we discuss also the detection prospects of the recycled pulsars in Sect. 5.

Fig. 8. Distribution in Galactic coordinates of the detected recycled pulsars in the simulations along with the observations.

4.5. Geometry

Examining the geometry of the pulsars detected in various wavelengths is interesting. Indeed, Fig. 9 shows the viewing angle, ζ, as a function of the inclination angle, χ for each detected pulsar in the simulation. Although these angles are generated between 0 and 180°, we folded them between 0 and 90°. Basically, we retrieve an area of detection slightly different for radio pulsars compared to Johnston et al. (2020) where they studied the normal pulsars, but almost the same area of detection for γ-ray pulsars. The radio pulsars are detected along a region much broader than along the ζ = χ diagonal as it was the case for normal pulsars, the pulsars in γ-ray surveys only are seen above the line of equation ζ = 90 − χ, and the pulsars detected in both radio and γ-ray surveys are the intersection of these two areas. The area of radio detection is broader for recycled pulsars compared to normal pulsars, because recycled pulsars have larger ρ in general, which broaden the area of detection |ζ − χ|≤ρ (or |ζ − (π − χ)| ≤ ρ). We remind that ∼71 % of the pulsars in γ-ray surveys only have a radio flux above 30 μJy and a geometry favorable for a radio detection. Thus many of them would still be seen with a follow-up observation in radio. One can notice that it seems particularly difficult to see a pulsar with a rotation axis and magnetic axis aligned, i.e. χ ≈ 0°. Indeed, for γ-ray detection, the condition |ζ − π/2|≤χ is needed, and below χ = 20°, it almost never happens. The rarity of low detected χ values appears to be due to the fact that, as shown in Yang & Li (2023), the effect of accretion tends to lead to a misalignment between the rotation axis and the magnetic axis, i.e. having χ increasing towards 90°. Therefore, after accretion, since the magnetic field is weaker, the torque n₃ caused by loss of energy via radiation of the pulsar has more difficulties to get χ to decrease. Furthermore, it is also plausible for those which get low χ that it took too much time and they became not luminous enough to be seen. These findings on the simulated population reflect well the observed values for ζ and χ. For instance Benli et al. (2021) constrained these angles for 10 radio-loud γ-ray MSPs, and found out that for both angles they are usually above 45°, and it is true that in our simulation the radio-loud γ-ray pulsars also verify this. Another example is the study of Guillemot & Tauris (2014), where they concluded that γ-ray energetic nearby MSPs are not seen because of their unfavorable viewing conditions, i.e. low ζ and we also see that high ζ are more favorable conditions of observations.

Fig. 9. Viewing angle (ζ) plotted against the inclination angle (χ) for the detected pulsars in the simulation.

4.6. Spider population

We may ask whether, at the end of the simulation, it is possible to identify spider pulsars within the recycled population. A key process governing the evolution of spider systems is the irradiation of the donor star by the pulsar wind as demonstrated in the works of Chen et al. (2013), Benvenuto et al. (2014), Misra et al. (2025). In particular, companions with masses between 0.1 and 1 M_⊙ can lose a significant fraction of their mass, causing the orbital period to first decrease and then increase again after a few Gyr (see Fig. 3 of Benvenuto et al. 2014, Figs. 1–4 of Chen et al. 2013, and Fig. 1 of Misra et al. 2025). We modelled this effect using the following equation:

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{c},\mathrm{evap}}=-f_{\mathrm{evap}}\frac{L_{\mathrm{sp}}}{2v_{\mathrm{c},\mathrm{esc}}^2}{\left(\frac{R_{\mathrm{c}}}{a}\right)}^2.\end{array}$$(33)

Indeed, the irradiation allows a mass loss Ṁ_c,evap from the companion, which depends on the escape velocity from the companion star surface v_c, esc, its radius R_c, the semi-major axis of the binary a, the pulsar spin-down luminosity L_sp and the efficiency parameter f_evap measuring how effectively the pulsar’s spin-down luminosity is converted into an evaporative wind. For our simulations we set f_evap = 0.05, as having greater values seem to make difficult the production of black widows systems (Chen et al. 2013). A typical diagnostic to distinguish redbacks (RBs) from black widows (BWs) is a plot of the binary orbital period (P_orb) versus the pulsar’s companion mass in M_⊙, as in the right panel of Fig. 10 with the observed systems (see also Fig. 1 of Blanchard et al. 2025). In Fig. 10, the left panel shows this diagram for the entire simulated population, while the right panel presents the same plot restricted to the detected systems in the simulation and observed systems. RBs are located in the upper right, and BWs appear in the left region of the plot. To compute P_orb in the simulation, we used Kepler’s third law, which gives

Fig. 10. Orbital period, Porb (in hr), as a function of companion mass (in M⊙). AGB stands for a companion that is on the asymptotic giant branch. Left panel: Entire simulated population. Right panel: Only systems with a detected pulsar by the pipeline in the simulation and observed systems. The companion mass reported is the median mass available in the ATNF catalogue for the observed systems. ‘Unknown (obs)’ means the pulsar has a null value in the ATNF catalog concerning its companion, ‘Sim’ stands for simulated, and ‘Obs’ stands for observed. The simulated systems are represented by squares and observed systems by circles. The different colours are for the different companions and are shown in the legend.

$$\begin{array}{c}\hfill P_{\mathrm{orb}}=\sqrt{\frac{4{\pi }^2{a}^3}{G(M_{\mathrm{c}}+M_{\mathrm{NS}})}},\end{array}$$(34)

where M_c is the companion’s mass, M_NS is the NS mass, a is the semi major axis of the binary and G is the gravitational constant. As a consequence, in a log-log plot of P_orb against M_c, a line with a −1/2 slope in the left panel for the BWs can be observed, as they also have a close orbit. Thus the dependence in a is diminished, which is contrary to the RBs that are in wider orbits. Observed RBs typically have main sequence companions. Observationally, there are 13 RB systems with a main-sequence companion, whereas our simulation produces nine such systems. We found slightly more pulsars in the RB region with COWD companions compared to the observations. Observed BWs have ultra-light companions, while most simulated systems in the BW region host either HeWDs or COWDs, suggesting that some of the companion of observed BWs could actually be HeWDs or COWDs, or we could discover such spiders. We also struggle to reproduce BWs systems with a P_orb > 2 h, whereas, on the left panel, we can see that a few systems with the characteristics of BWs with P_orb > 2 h are formed but not detected. In addition, there are currently around 32 confirmed RBs and 50 confirmed BWs (Koljonen & Linares 2025), whereas our simulation produces roughly three times fewer spider-like pulsars. These discrepancies may indicate that a refinement is necessary for the initial conditions concerning the semi-major axis a which determine the initial orbital period, or combining a different set of parameters for the efficiency parameter of irradiation f_evap and the accretion efficiency f_MT or even other parameters describing the binary evolution. It may also hints that these systems might often form through alternative channels, as proposed by Metzler & Wadiasingh (2025). These authors suggest that the ultra-light companions of spiders could form in high-metallicity accretion discs over timescales exceeding > 1 Gyr (Patruno & Kama 2017), through dynamical capture or exchange interactions in clusters (Sigurdsson et al. 2003), or, especially for BWs, by the accretion of an ultra-compact X-ray binary companion’s outer layers onto the NS (Bailes et al. 2011). The study by Misra et al. (2025), which uses the detailed stellar evolution code MESA to characterise the spider population, shows that an accretion efficiency of 0.7 works better to reproduce the observed systems. We did not find any major differences in terms of detections between a simulation with f_MT = 0.7 or also by changing f_evap = 0.5. In summary, the simulation successfully produces two distinct groups reminiscent of spiders, which is encouraging. However, it underestimates the total number of observable spider systems. Therefore, including the additional formation channels mentioned above and fine-tuning the parameters related to binary evolution would likely lead to a better agreement with observed populations.

4.7. AIC scenario

Another formation channel for pulsars is the accretion-induced collapse (AIC) of WDs (Wang & Gong 2023). We performed a simulation starting with one million binary systems composed of an ONeWD and a main sequence companion, with the latter at a random evolutionary stage (0–100% of its main sequence lifetime), as was done for the pulsars in this study. We considered only ONeWDs, as their higher mass and composition make them ideal candidates for electron-capture supernovae. Consequently, most of the evolutionary channels leading to AIC involve an ONeWD (Wang & Liu 2020). Additional simulations were conducted with an initial population consisting of one third of each WD type, and only the ONeWDs were found to undergo collapse into NSs. The initial mass distributions for the WDs were taken from Tremblay et al. (2016). Figure 11 shows the resulting P–Ṗ diagram: the left panel displays the entire population of NSs formed via AIC, while the right panel presents only those that are detectable (either in radio or γ-rays) based on our detection pipeline, alongside the observed pulsar population. Out of the one million initial systems, ∼9000 WDs evolved into NSs (∼0.9%). Within the total population, we distinguish two main groups: one consists of NSs that accreted little mass and remain relatively young, behaving like “normal” pulsars; the other group lies in the region of mildly recycled pulsars (those with 30 ≤ P ≤ 1000 ms and Ṗ < 10⁻¹⁶ s/s), having accreted significantly more mass and consequently spun up to shorter periods, although it was not sufficient for them to become fully recycled as MSPs. After applying the detection pipeline, only ∼600 pulsars remain. Most of these belong to the normal pulsar population. A few objects (∼70) are found in the mildly recycled region of the P–Ṗ diagram, but they exhibit relatively high spin periods, with none below 20 ms. This suggests that a non-negligible fraction of recycled pulsars with relatively high spin periods may originate from the AIC of ONeWDs (depending on the birth rates of these objects). As shown in Fig. 3, these long-period/mildly recycled pulsars are not fully accounted for by the formation channel involving an NS formed via core-collapse supernova that later accretes mass from a companion. This discrepancy explains the imperfect agreement of our KS test for the radio population. The existence of this alternative formation pathway could help alleviate the so-called birth rate problem (Gençali & Ertan 2024), in which the Galactic core-collapse supernova rate appears to be lower than the combined estimated birth rates of NSs. In summary, this informs that among the whole ONeWDs simulated, ∼0.9% collapse into NSs, ∼6.9% are detectable and then ∼12.2% are in the area of mildly recycled pulsars considered in this work. Therefore, ∼7.6 × 10⁻³% of ONeWD (that were in a binary system) seem to contribute to the population of mildly recycled pulsars. The contribution of ONeWDs to the pulsar population could be explored in depth in future work, as it would require knowledge of the birth rate of ONeWDs, a quantity that remains highly uncertain.

Fig. 11. Diagram of P−Ṗ of pulsars born after the AIC of WDs. Left panel: Entire simulated population. Right panel: Only systems with a detected pulsar by the pipeline along with the observations.

It was suggested in Tauris et al. (2013) that the observed pulsars located in region II of the Corbet diagram (Fig. 12) may originate from AIC. In Fig. 12, the recycled pulsars originating from an AIC channel are shown as red circles with black outlines, whereas the other points correspond to recycled pulsars formed through core-collapse supernovae, paired either with a HeWD or a COWD companion. In our simulation, all pulsars formed through AIC indeed end up in region II. However, we also produce pulsars in region II that do not have an AIC origin. Therefore, while identifying an observed system as AIC-formed would imply that it should lie in region II, the reverse is not true: a system found in region II does not necessarily result from AIC. That said, if a system has P_orb < 1 day and a P > 10 ms, an AIC origin appears significantly more likely.

Fig. 12. Spin period (P) plotted against the orbital period (Porb) of the simulated pulsars. The red circles with a black outline are the AIC WD that became recycled pulsars, while the other points represents recycled pulsars formed via core-collapse supernova with either a HeWD or COWD companion. The four regions, labelled I, II, III, and IV, represent the same cut as in Tauris et al. (2013) to explain their differences as an evolutionary point of view.

5. Discussion

5.1. Model evaluation and caveats

The model we used seems to reproduce the recycled pulsars in the radio and γ-ray wavelengths. However, it is important to address several points. A crucial aspect of the model is the mass transfer. Changing the accretion efficiency, f_MT, was not explored in this work, and we decided to keep the value proposed by SEVN and used in other similar works (Bouffanais et al. 2021; Iorio et al. 2023). Nevertheless, changing this value would clearly impose a change on some parameters, such as τ_d and the birth rate.

Furthermore, SEVN does not allow spin-up during the CE phase. Therefore, this is a channel of formation of recycled pulsars from which we do not evaluate the contribution. However, Osłowski et al. (2011) estimated that the accreted matter during the CE phase would be too chaotic to produce spin-up and it could only decay the magnetic field. We also do not take into account the possible accretion from winds. Kiel et al. (2008) estimated that the spin-up caused by wind would also be too chaotic to produce any significant spin-up.

Still regarding the CE phase, major uncertainties remain about its treatment in binary evolution (Tauris & van den Heuvel 2023), particularly for the formation of double neutron star (DNS) or neutron star–black–hole (NS+BH) systems. When the donor star fills its Roche lobe, a CE phase is triggered with the NS, leading to a strong reduction of the orbital separation before the donor undergoes its supernova explosion. Moreover, Kruckow et al. (2016) showed that an in-spiraling NS can in principle eject the envelope of its companion, but that predicting the final post-CE separation remains highly uncertain. As long as dedicated studies do not clarify the CE evolution for such systems, accurately reproducing the number of observed DNS and NS+BH will remain challenging. In our simulations, we form about ∼2 − 5 DNS or NS+BH systems, while observations report 12 DNS systems. Despite these uncertainties, our models are still able to produce a fraction of the observed population. However, the small number of systems involved implies that the statistical significance remains limited.

The parameter ϵ characterising the ellipticity of the NS could also be an interesting parameter to vary in order to explore the possibility of the formation of sub-MSPs (Haensel & Zdunik 1989). In this work, we only tried to draw this parameter from a uniform distribution between the minimum and maximum value of ϵ hypothesised by Giliberti & Cambiotti (2022). The works of Bildsten (1998), Chakrabarty et al. (2003) indeed suggested that GW radiation allows the spin period P to remain above the millisecond threshold. Therefore, giving the same ϵ to every NS in the simulation and trying to find the value to obtain sub-MSPs could be quite interesting. Everything that is mentioned above could be explored in future works.

Part of the angular momentum of the binary is lost through magnetic braking, and one as to know that we used the formula proposed by Rappaport et al. (1983) for this phenomenon as it is implemented by default in SEVN. However, recently Chen et al. (2021) showed that using the prescription of Van et al. (2019) of magnetic braking could, among other things, influence a lot the number of NS + WD systems. Thus, in the future, it would also be interesting to evaluate the influence on the population while modifying the magnetic braking model of the binary.

For the simulation involving ONeWDs and main sequence stars, we adopted the same initial distributions for eccentricity (e) and mass ratio (q) as used for massive binaries, as described in Sect. 2. We also explored alternative power-law indices for the q distribution, varying it between −0.5 and 0.5 instead of the default −0.1, but found no significant impact on the simulation outcomes. Similarly, due to the lack of reliable constraints on the birth eccentricity distribution for low-mass binaries, we tested different power-law indices for e as well, and again observed negligible differences in the results.

Let us mention several geometrical or physical limitations concerning the striped wind model, which may affect our results. For instance, as suggested by Kalapotharakos et al. (2022), a more accurate description of the γ-ray luminosity might involve an additional parameter, such as the cut-off energy, and a different dependence on B and P, or equivalently on P and Ṗ, than the one assumed in (29). Moreover, plasma recollimation effects within the striped wind could reduce the γ-ray beam opening angle compared to the adopted value, leading to a smaller number of detectable pulsars.

5.2. Comparison with other works

Recently, Tabassum & Lorimer (2025) investigated the populations of both normal pulsars and MSPs in the radio and γ-ray bands. Their study focused in particular on testing various γ-ray luminosity models, leading them to rule out three of them–namely, the pair-starved polar cap (PSPC), the slot gap two-pole caustic (TPC), and the outer gap (OG) models–because they overproduced detectable pulsars. However, their overall modelling framework differs significantly from ours: they do not simulate the full accretion process, they do not include magnetic field decay, and they assume an exponential decay of the inclination angle χ (referred to as α in their work). Moreover, when they adopt the same γ-ray luminosity prescription as in this study, they do not incorporate the striped wind geometry. As a result, their model predicts a much larger number of detectable γ-ray sources, ranging from 350 up to 650. In contrast, our more comprehensive approach, which includes magnetic field decay, accretion history, and the striped wind geometry, yields a tighter prediction, with between 290 and 330 detectable γ-ray recycled pulsars, reflecting stochastic variations across simulation realisations. This results in fewer than ∼220 unidentified sources in the 4FGL catalog, significantly narrowing the allowed range compared to the upper limit of ∼520 unidentified sources predicted by Tabassum & Lorimer (2025).

Gonthier et al. (2018) estimated that approximately 11 000 MSPs would be required within the inner Galactic centre (defined as the 7° ×7° region of interest, corresponding to a radius of about 0.5 kpc) to account for the GeV excess. Their scenario relies on the contribution of GCs that have migrated into the Galactic bulge through dynamical friction (Gnedin et al. 2014). In contrast, our study considers only pulsars formed in the Galactic spiral arms. These tend to be more luminous than those assumed by Gonthier et al. (2018). Nevertheless, our results indicate that only about 190 MSPs born in the spiral arms contribute to roughly 5% of the GeV excess. Therefore, if the GeV excess is entirely due to pulsars, it is still entirely plausible that a much larger population of MSPs originating from GCs resides in the Galactic centre as estimated by Gonthier et al. (2018) or they were born in the bulge, but they are beyond current detectability.

The studies of Kiel et al. (2008), Kiel & Hurley (2009) were among the first to attempt reproducing the pulsar population by accounting for the full evolutionary path of NSs in binary systems. However, several key differences exist between their approach and ours. First, their magnetic field decay prescription is exponential; thus, if they had used the same decay timescales τ_d as those identified in our work (0.05, 1.4, or 2.3 Myr), the magnetic fields of their NSs would decay unrealistically fast. As a result, they require significantly higher values of τ_d to reproduce the observed population. Secondly, the equation they adopt for the evolution of the angular velocity Ω is much simpler than the one used in this study, which incorporates several additional physical effects. Furthermore, while our initial distributions for P₀ and B₀ are informed by observational constraints and recent population synthesis efforts (Igoshev et al. 2022; Sautron et al. 2024), their work relies on uniform distributions for these parameters, not grounded in observational data or detailed modelling. Finally, their comparison with observations is limited, as their population synthesis framework does not include a detection model to account for selection effects.

As mentioned in Sect. 2, we obtained matching P−Ṗ diagrams from the observations and simulations for a number of simulated binaries between 4.6 × 10⁵ and 5.7 × 10⁵, corresponding to a birth rate between 3.3 − 4.1 × 10⁻⁶ yr⁻¹. Studies by Story et al. (2007) and Gonthier et al. (2018) found the birth rate of MSPs to be in the range of 3 − 6.5 × 10⁻⁶ yr⁻¹. The birth rate found in our study is consistent with the findings of Story et al. (2007) and Gonthier et al. (2018). We constrain the birth rate of recycled pulsars even more with this work.

5.3. Detection prospects

We wanted to assess the detectability of new recycled pulsars sources with this model. To do so, we considered the survey parameters of SKA. See Appendix D to see the parameters and in order to make predictions for radio observations. Moreover, concerning the γ-ray wavelengths, we took into account an instrument that would be ten times more sensitive than Fermi/LAT (i.e. we adopted a detection threshold of F_lim = 1 × 10⁻¹⁶ W.m⁻², as the average flux limit of ∼1 × 10⁻¹⁵ W.m⁻² was derived from the Fermi/LAT all-sky sensitivity map). These considerations allowed us to multiply the number of radio detections by approximately three and the number of γ-ray detections by about five. This improvement in sensitivity also allowed us to detect around ten γ-ray pulsars in the Galactic centre, indicating that an even higher sensitivity would be necessary to uncover additional pulsars in the Galactic centre.

6. Conclusions

This study is the first ever to evolve the population of pulsars in binary systems in a complete way by taking into account the binary processes. By doing so, we also took into account the observational biases from radio and γ-ray surveys in order to analyse the recycled pulsar population.

We showed that our model reproduces most of the recycled pulsar population in the radio and γ-ray wavelengths and that most of the mildly recycled pulsars (those with 30 ≤ P ≤ 1000 ms and Ṗ < 10⁻¹⁶ s/s) seem to originate from the AIC of WDs. In order to reproduce the population, the characteristic B-field decay timescale, τ_d, was chosen to be 0.05, 1.4, or 2.3 Myr, with most of the population evolving with τ_d = 0.05 Myr. Moreover, with this model we estimate a birth rate for recycled pulsars between (3.3–4.1) × 10⁻⁶ yr⁻¹. The striped wind model allowed us to reproduce the gamma-ray peak separation observed for γ-ray recycled pulsars very well. Most of the pulsars accrete for 10^8 − 9 yr no matter the type or initial mass of the companion. Since the timescale for alignment between the rotation axis of the NS and the normal of the orbital plane is often short, we obtained ∼80% of the population that has these two axes aligned (Lorange et al. 2026). Moreover, the mean mass obtained for this population is about 1.8 M_⊙, with some very high mass pulsars, up to 2.7 M_⊙. The mean mass is in agreement with the recent measurements of high masses of pulsars in spider systems (e.g. Linares 2020). Furthermore, in the simulation we were able to recover pulsars that have the characteristics of spiders but not the whole observed population, indicating that other channels of formation (in addition to the one considered in this work) must be at work to obtain spiders. Moreover, to obtain a more similar population of spiders, as suggested by Misra et al. (2025), it is likely necessary for a fraction of the population to have a higher accretion efficiency. The positions of the detected recycled pulsars in the simulation are generally well reproduced, although the simulation predicts more pulsars between latitudes −15° and 15° than observed. We estimate that MSPs formed in the Galactic field contribute only to ∼5% of the γ-ray galactic center excess, corresponding to roughly 190 MSPs in the Galactic centre. This result suggests that if the γ-ray galactic center excess originates only from pulsars, the majority of MSPs in the Galactic centre likely come from GCs that migrated inwards through dynamical friction or they were born in the bulge. Either way, we are not able to detect them today (Gnedin et al. 2014; Gonthier et al. 2018). Finally, most of the detected pulsars seem to have an inclination angle, χ, greater than 10°, and it also seems easier to observe pulsars that have a viewing angle, ζ, greater than 45°, which has been supported by other studies (Benli et al. 2021; Guillemot & Tauris 2014). We were also able to make some predictions: We expect approximately three times more radio recycled pulsars thanks to SKA and about five times more detections of γ-ray recycled pulsars with an instrument that has a sensitivity increased by one order of magnitude compared to Fermi/LAT.

Other interesting topics to explore include adding another channel of formation for recycled pulsars; studying the contribution of AIC ONeWD to the population of pulsars in depth; adding the possible ablation of the companion by the pulsar to study the spiders; and evaluating the detectability of continuous GWs from the recycled population, as was already done for the normal population by Cieślar et al. (2021). In addition, more constraints could be found about the numerous parameters of the model while focusing on a restricted number of parameters.

Acknowledgments

This work has been supported by the French Research Agency grant ANR-20-CE31-0010. We would like to thank Francesca Calore, Joanna Berteaud and Maïca Clavel for insightful discussions. D.M. acknowledges the support of the Department of Atomic Energy, Government of India, under project No. 12-R&D-TFR-5.02-0700. We also thank the anonymous referee for the helpful comments and suggestions, which helped improve the paper.

References

Appendix A: Ratio of the pulsar’s radiation loss torque and the gravitational wave torque

In the set of differential equations that we use (Eqs. (9a) to (9c)), we take into account a torque, n_GW, caused by GWs. However, this torque depends on Ω⁵. Therefore compared to the torque due to the pulsar’s radiation loss, n₃ (which depends to Ω³), it can often be neglected. If we compute the ratio of both torques, we obtain

$$\begin{array}{c}\hfill \frac{n_{\mathrm{GW}}}{n_3}=\frac{G{I}^2{\mathrm{\Omega }}^2{ϵ}^2{\mu }_0}{5{\pi }^3{c}^2{B}^2{R}^6}.\end{array}$$(A.1)

To choose when to consider the torque caused by GWs, we look into when the ratio of Eq. (A.1) would be greater than 0.05. If we take B = 7 × 10⁴ T and ϵ = 10⁻⁷, it gives Ω_lim ∼ 1100 rad.s⁻¹, while actually most of the time ϵ can be lower than this, or the magnetic field a bit higher, resulting in higher Ω_lim. Thus, to be very cautious and thorough, we even adopted Ω_lim = 960 rad.s⁻¹. Therefore in this work, we do not neglect the GWs torque for Ω ≥ Ω_lim.

Appendix B: Influence of the ISM on the radio pulse profile

Dispersion and scattering caused by the ISM lead to a broader radio pulse profile. In order to determine the influence of the ISM through τ_DM, the formula of Bates et al. (2014) can be used,

$$\begin{array}{c}\hfill {\tau }_{\mathit{DM}}=\frac{e^2\mathrm{\Delta }{f}_{\mathrm{ch}}DM}{4{\pi }^2{ϵ}_0{m}_{e}c{f}^3}=8.3\times {10}^{15}s\frac{\mathrm{\Delta }{f}_{\mathrm{ch}}}{f^3}DM.\end{array}$$(B.1)

Here, e is the electron charge, m_e its mass, Δf_ch the width of frequency of the instrument channel in Hz, f the observing frequency in Hz, and DM the dispersion measure in pc/cm³. Furthermore, in order to determine the influence of scattering by an heterogeneous and turbulent ISM through τ_scat, the empirical fit relationship from Krishnakumar et al. (2015) is used

$$\begin{array}{c}\hfill {\tau }_{\mathrm{scat}}=3.6\times {10}^{-9}D{M}^{2.2}(1+1.94\times {10}^{-3}D{M}^2).\end{array}$$(B.2)

Both τ_DM and τ_scat are given in units of second (s) and depend on the dispersion measure DM (in units of pc/cm³), which is found for each pulsar by running the code of Yao et al. (2017) converting the distance of a pulsar into the dispersion measure thanks to their state-of-the art model of the Galactic electron density distribution. See Table 1 for the parameters used for the two surveys considered, i.e. FAST GPPS and PMPS.

Appendix C: Anisotropy of the striped wind model

The link between the γ-ray flux and the total γ-ray luminosity is not a straightforward relation to deduce. The problem resides in the highly anisotropic emission pattern of the pulsar striped wind. This anisotropy is usually embedded in a so-called beaming factor, f_Ω, describing the flux pattern depending on the obliquity χ and line of sight inclination ζ. Its definition is found in Watters et al. (2009), Pétri (2011, 2026) and shown on a log scale in fig. C.1 for the force-free dipole magnetosphere model. For high inclinations (ζ), it is about 0.2, whereas for low inclinations( ζ), it is very high, well above ten.

Fig. C.1. Flux anisotropy factor (fΩ; in log) for the striped wind model versus obliquity (χ) and inclination of Earth’s line of sight (ζ).

Appendix D: SKA survey parameters

Table D.1 displays the SKA survey parameters in the mid-frequency ranges we used to make the predictions in Sect. 5.3. SKA is described in Carilli & Rawlings (2004), Kramer & Stappers (2015), Grainge et al. (2017). SKA-1-Mid represents the estimate of the initially planned SKA operation.

All Tables

All Figures

	Fig. 1. Diagram of P−Ṗ of the observed pulsars considered in this work along with the death line (solid green line) and the death valley (shaded green area).
In the text

	Fig. 2. Evolution in the P−Ṗ diagram of a pulsar in the simulation that becomes an MSP.
In the text

Fig. 3. Diagram of P−Ṗ for observed and simulated pulsars in different wavelengths. Left panel: P−Ṗ diagram for the detected pulsars in radio surveys only in the simulations in red along with the observations in radio for the FAST GPPS or PMPS survey in blue. Middle panel: P−Ṗ diagram for the detected pulsars in γ-ray surveys only in the simulations in green, along with the observed pulsars in γ-ray from Fermi surveys only in blue. Right panel: P−Ṗ diagram for the simulated pulsars detected both in radio and γ-ray surveys in purple, along with observations for the FAST GPPS or PMPS survey simultaneously with the γ-ray observations from Fermi in blue. For all the panels, the death line is represented as a green solid line, and the death valley is the shaded green area.

In the text

	Fig. 4. Distribution of the γ-ray peak separation of the simulated population along with the observations from the 3PC catalogue.
In the text

	Fig. 5. Width of the radio profile plotted against the spin period, P, of the radio-detected pulsars in the simulation along with the observations.
In the text

Fig. 6. Companion outcomes in the pulsar population. Left panel: Distribution of the type of the pulsars companions. HeWD stands for helium white dwarf, ‘isolated’ means the pulsar does not have a companion anymore, COWD stands for carbon-oxygen white dwarf, NotARemnant means the companion is a non-degenerate star, NS CCSN indicates an NS born after a core-collapse supernova, NS ECSN refers to an NS born after an electron capture supernova, and ONeWD stands for oxygen-neon white dwarf. Right panel: Initial mass of the pulsars’ companions plotted against the duration of the accretion phase, with the final type of the companion displayed. MS stands for main sequence star, and NS stands for NS.

In the text

	Fig. 7. Histogram of the final mass of the recycled pulsars in the simulation. The vertical dashed black line represents the theoretical maximum mass reachable by an NS from the study of Ruiz et al. (2018).
In the text

	Fig. 8. Distribution in Galactic coordinates of the detected recycled pulsars in the simulations along with the observations.
In the text

	Fig. 9. Viewing angle (ζ) plotted against the inclination angle (χ) for the detected pulsars in the simulation.
In the text

Fig. 10. Orbital period, Porb (in hr), as a function of companion mass (in M⊙). AGB stands for a companion that is on the asymptotic giant branch. Left panel: Entire simulated population. Right panel: Only systems with a detected pulsar by the pipeline in the simulation and observed systems. The companion mass reported is the median mass available in the ATNF catalogue for the observed systems. ‘Unknown (obs)’ means the pulsar has a null value in the ATNF catalog concerning its companion, ‘Sim’ stands for simulated, and ‘Obs’ stands for observed. The simulated systems are represented by squares and observed systems by circles. The different colours are for the different companions and are shown in the legend.

In the text

	Fig. 11. Diagram of P−Ṗ of pulsars born after the AIC of WDs. Left panel: Entire simulated population. Right panel: Only systems with a detected pulsar by the pipeline along with the observations.
In the text

Fig. 12. Spin period (P) plotted against the orbital period (Porb) of the simulated pulsars. The red circles with a black outline are the AIC WD that became recycled pulsars, while the other points represents recycled pulsars formed via core-collapse supernova with either a HeWD or COWD companion. The four regions, labelled I, II, III, and IV, represent the same cut as in Tauris et al. (2013) to explain their differences as an evolutionary point of view.

In the text

	Fig. C.1. Flux anisotropy factor (fΩ; in log) for the striped wind model versus obliquity (χ) and inclination of Earth’s line of sight (ζ).
In the text