On the Formation of Double Neutron Stars in the Milky Way: Influence of Key Parameters

The detection of gravitational wave events has stimulated theoretical modeling of the formation and evolution of double compact objects (DCOs). However, even for the most studied isolated binary evolution channel, there exist large uncertainties in the input parameters and treatments of the binary evolution process. So far, double neutron stars (DNSs) are the only DCOs for which direct observations are available through traditional electromagnetic astronomy. In this work, we adopt a population synthesis method to investigate the formation and evolution of Galactic DNSs. We construct 324 models for the formation of Galactic DNSs, taking into account various possible combinations of critical input parameters and processes such as mass transfer efficiency, supernova type, common envelope efficiency, neutron star kick velocity, and pulsar selection effect. We employ Bayesian analysis to evaluate the adopted models by comparing with observations. We also compare the expected DNS merger rate in the galaxy with that inferred from the known Galactic population of pulsar-neutron star systems. Based on these analyses we derive the favorable range of the aforementioned key parameters.


INTRODUCTION
Since the discovery of the first chirp of gravitational waves (GWs) from a binary black hole (BH) merger GW 150914 by the Laser Interferometer Gravitational-wave Observatory (LIGO; Aasi et al. 2015;Abbott et al. 2016), there have been 90 compact binary merger events detected (The LIGO Scientific Collaboration et al. 2021).The first double neutron star (DNS) merger GW event GW170817 (Abbott et al. 2017a) was found to be accompanied by a short gamma-ray burst (SGRB, Abbott et al. 2017b) and a kilonova (Abbott et al. 2017c), providing direct evidence for the connection between these phenomena.
The detection of GW events has also stimulated theoretical modelling of the formation and evolution of double compact objects (DCOs).The proposed formation channels for merging DCOs include isolated binary evolution, chemically homogeneous evolution, dynamical interactions in dense environments, hierarchical triple evolution, binary evolution in active galaxy nucleus disks, etc (see Barack et al. 2019;Mandel & Broekgaarden 2022;Mapelli et al. 2022, for recent reviews).For the first time, it is possible to quantitatively compare the merging rate predictions against observations.However, ever for the most studied isolated binary evolution channel, there exist big uncertainties in the input parameters and in the treatments of the binary evolution processes.For example, a large fraction of merging DCOs have experienced at least one common envelope (CE) evolution phase (Ivanova et al. 2013, for a review).Depending on the CE efficiency α and the binding energy parameter λ of the donor's envelope, this channel can produce a large range of the merging rates.Since they can be broadly consistent with the merger rate inferred from GW observations, it is difficult to put stringent constraints on the above mentioned parameters.In particular, some of the models still adopt a constant λ, which seems physically unrealistic.
So far, DNSs are the only DCOs for which direct observations are available through traditional electromagnetic astronomy.A viable model predicting the properties of merging DNSs detectable by GW observatories should also be simultaneously able to account for the Galactic DNS population, although the known 21 DNS pulsars from radio observations represent a very small fraction of all DNSs in the Galaxy.Phinney (1991) and Narayan et al. (1991) tried to infer the DNS merger rates from binary pulsar observations.Since then, there has been remarkable progress in both observation and theory including the growing DNS data and increasingly accurate treatments of the selection effects.
In this work, we construct a large number of models for the formation of Galactic DNSs, taking into account various possible combinations of the critical input parameters and processes like mass transfer efficiency, supernova (SN) type, CE efficiency, NS kick velocity, and selection effect of pulsars.Through both qualitative and quantitative comparison with the Galactic DNSs, we aim to obtain useful constraints on the parameters and understanding of the binary evolution processes.Here we focus on isolated binary evolution, since most DNSs are in the Galactic field (see below) and dynamical formation of DNSs in star cluster environments does not contribute significantly to the overall merger rate (Phinney 1991;Belczynski et al. 2018;Fong et al. 2019;Ye et al. 2019).
Since the first DNS system PSR B1913+16 was discovered by Hulse & Taylor (1975), there are now 21 DNSs known in the Galaxy (see Table 1, part of data are from the ATNF pulsar catalog1 , Manchester et al. 2005).Among them, 18 DNSs are in the Galactic field and 3 are in globular clusters (GCs).The Galactic field DNSs are thought to be predominantly formed through isolated binary evolution (Tauris et al. 2017).
Most of the pulsars in DNSs are the first-born NS, and are observed as (mildly) recycled pulsars.They have accreted more or less mass from the progenitor star of the secondly born NS.Because of mass accretion-induced field decay (Bhattacharya 2002), the recycled pulsars have a weaker magnetic field than normal pulsars, and a longer observable duration in radio (see Fig. 2 of Chattopadhyay et al. 2020).Therefore, we limit our study of the radio characteristic of DNSs to the first-born pulsars.
In this paper, we combine binary population synthesis (BSE) with pulsar evolution to explore the formation of Galactic DNSs, their radio characteristics, and expected merger rates.The aim of this work is to evaluate the influence of the key parameters by comparing theory with observation.The remainder of this paper is organized as follows.We introduce isolated binary evolution, the BSE code and the pulsar evolution model in Section 2. We demonstrate the calculated results in Section 3, and explore the influence of the radio selection effect in Section 4. We present our discussion and conclusions in Section 4.

Binary Population Synthesis
We first briefly introduce the two dominant channels of DNS formation by isolated binary evolution (Bhattacharya & van den Heuvel 1991;Tauris & van den Heuvel 2006;Vigna-Gómez et al. 2018).In the first channel (Channel I), the evolution begins from a primordial binary with both components more massive than 8M ⊙ , and the mass ratio (q = primary mass/secondary mass) considerably larger than unity.The primary star evolves more rapidly to overflow its Roche-lobe (RL) and transfers mass to the secondary.After its envelope is exhausted, the primary evolves to be a helium star, and finally a NS with a SN explosion.The subsequent evolution of the secondary star leads to a CE phase, in which the NS spirals into the envelope of the secondary star.If the NS is able to expel the envelope, the outcome is a NS-helium star binary in a compact orbit.A DNS finally forms if the binary is not disrupted during the second SN.In the second channel (Channel II), the initial mass ratio q ∼ 1, and the binary experiences a double-core CE phase.This may leave a double helium star binary in a compact orbit.The DNS then evolves from the double helium star binary (Brown 1995;Dewi et al. 2006;Hwang et al. 2015;Andrews et al. 2015).
For the isolated binary evolution, there are some critical physical processes that have not been well understood, such as the stability of mass transfer (MT), the fate of CE evolution, the SN mechanisms and the natal kicks to the NS (e.g., Andrews et al. 2015;Vigna-Gómez et al. 2018;Shao & Li 2018a;Chu et al. 2022).They will be discussed in some detail below.Besides them, before the second SN, there may exist a Case BB MT phase during which the binary comprises a NS and a naked helium star companion (Habets 1986).The Case BB phase also plays an important role in the formation of DNSs (Vigna-Gómez et al. 2018).During this phase, the first-born NS undergoes mass accretion and is accelerated to short spin period (∼ 10 − 100 ms, depending on the accretion rate and the amount of accreted mass).This may result in an ultra-stripped helium star before the SN, which will finally evolve to be an ONeMg WD or a NS depending on the initial orbital separation and the mass of the helium star (Tauris et al. 2013(Tauris et al. , 2015(Tauris et al. , 2017)).
After the birth of the DNS, its orbit keeps shrinking by GW radiation (Peters 1964;Hulse & Taylor 1975), finally producing a merger event that may be associated with a SGRB (depending on the assumptions about the DNS masses and other factors, Abbott et al. 2017a;Mösta et al. 2020;Mandel & Broekgaarden 2022).The merger time (τ gwr ) is determined by the orbital separation, eccentricity and the NS masses.Andrews & Mandel (2019) accordingly split the Galactic field DNSs into three categories: subpopulation (i) have small eccentricities and τ gwr less than Hubble time; subpopulation (ii) contain widely separated systems and will not merge within Hubble time; and subpopulation (iii) have large eccentricities and τ gwr less than Hubble time (see also Fig. 1).
We use the BSE code originally developed by Hurley et al. (2002) to calculate the Galactic DNS distribution.We adopt a star formation rate of 3 M ⊙ yr −1 , and follow the binary evolution over 12 Gyr.The primordial binaries consist of a primary star with mass M 1 and a secondary star with mass M 2 in a circular orbit with separation a (or orbital period P orb ).For the distributions of M 1 , M 2 and a, we set n χ grid points of the parameter χ logarithmically spaced with δlnχ = 1 n χ − 1 (lnχ max − lnχ min ). (1) If a specific binary i evolves to a DNS, then the binary contributes to the DNS population with a rate where f bin is the binary fraction (we assume all-stars are initially in binaries, i.e., f bin = 1.0),SF R is the star formation rate, and M ⋆ = 0.5M ⊙ is the average mass for all stars; Φ, ϕ, and Ψ are the distribution functions.
We assume the distribution of the primary masses to follow the initial mass function ξ(M 1 ) of Kroupa et al. (1993) in the range of 5 − 40M ⊙ .Thus, Φ(lnM 1 ) = M 1 ξ(M 1 ). (3) For the secondary masses, we assume that they follow a flat distribution between 0 and M 1 (Kobulnicky & Fryer 2007), then, The initial orbital separations a are taken to be uniformly distributed in the range of 3 − 10 4 R ⊙ in the logarithm (Abt 1983), so Ψ(ln a) = k = 0.12328. (5) Besides the standard treatments of stellar and binary evolution in the BSE code, we modify and update some key processes in the formation of DNSs.They are described as follows.

Mass loss, mass transfer and CE evolution
We adopt the fitting formula of Nieuwenhuijzen & de Jager (1990) to model stellar wind mass loss.For hot OB stars (T eff > 11000 K) and stripped helium stars, we replace it with the simulated relations by Vink et al. (2001) and Vink (2017), respectively.
When the primary star evolves to fill its RL, the first MT begins.The accretion efficiency play a vital role in determining the fate of the binary evolution.If the MT is sufficiently rapid, the accreting secondary star will be out of thermal equilibrium, and the responding expansion may lead the secondary star to fill its own RL, resulting in a contact binary (Nelson & Eggleton 2001).In this situation, if the primary star is a main sequence star, the binary will merge; otherwise, it will enter a CE phase (de Mink et al. 2013).
Here we adopt three models for the accretion efficiency in the primordial binaries.Model I assumes that the mass accretion rate Ṁ2 of the secondary star is modulated by its rotation rate: Ṁ2 = Ṁ1 × (1 − Ω2 Ωcr ), where Ṁ1 , Ω 2 , and Ω cr are the MT rate, the angular velocity of the secondary and its critical value (break-up rate), respectively.The rotation of the star is set to be rigid, and its variation is affected by both accretion and tidal synchronization2 .
Model II assumes that half of the transferred mass is accreted by the secondary.Model III assumes that the transferred mass is always accumulated by the secondary unless its thermal timescale becomes significantly shorter than the mass transfer timescale.So Ṁ2 is set to be the MT rate multiplied by a factor of min(10 τ Ṁ τKH,2 ,1) where τ Ṁ is the MT timescale and τ KH,2 is the secondary's thermal timescale (Hurley et al. 2002).Rapid mass accretion may push the accretor beyond its thermal equilibrium, causing it to expand and become overluminous.According to the results of Shao & Li (2014), τ KH,2 is usually less than τ Ṁ , which means that the MT is almost conservative.
In all the three models, the excess material is assumed to be ejected out of the binary in the form of isotropic winds, taking away the accretor's specific angular momentum.Thus, the accretion efficiency affects not only the orbital evolution of the binary system, but also the stability of MT (Soberman et al. 1997).According to the calculation by Shao & Li (2014), the maxium q cr for stable MT are ∼ 6, 2.5, and 2.2 for models I, II and III, respectively.We adopt their results (Fig. 1) of q cr to deal with the MT processes in the primordial binaries.
In the case of NSs accreting from a main-sequence/giant donor star, we use the Eddington limit to constrain the accretion rate, and adopt the maximal mass ratio q cr = 3.5 for stable MT following previous works (Tauris et al. 2000;Podsiadlowski et al. 2002;Shao & Li 2012;Deng et al. 2020).For NS-Helium star binaries, we assume that the MT becomes dynamically unstable when the orbital period is less than 0.06 day (Tauris et al. 2015).
To deal with CE evolution, we compare the orbital energy of the binary and the binding energy of the donor's envelope, to evaluate whether the orbital energy is sufficient to expel the envelope (Webbink 1984), where α CE is the CE efficiency, M 1c is the mass of the primary's core, a i and a f are the orbital separations of the binary before and after the CE stage respectively, G is the gravitational constant, r L is the RL radius of the core, and λ is the binding energy parameter of the primary's envelope.We let α CE = 0.1, 0.5, 1, 2, 5 and 10 to evaluate its possible influence.In the case with α CE > 1 it is assumed that there are extra energy sources (such as recombination energy) used to expel the stellar envelope (values significantly larger than 1 might not reflect real physical process but for illustrative purposes).The binding energy parameter λ depends on the stellar evolutionary stage, and we adopt the calculated results of Xu & Li (2010a,b) and Wang et al. (2016).

SN models and natal kicks
We use the rapid and delayed explosion models proposed by Fryer et al. (2012) and the stochastic explosion model proposed by Mandel & Müller (2020) to process the products of core-collapse supernovae (CCSN).In the rapid and delayed SN models, the CO core masses at the SN determines the remnant masses, which are further increased by accreting from the subsequent fallback material.Stars with inital mass ∼ 9 M ⊙ are thought to produce electroncapture supernovae (ECSN), and the helium core mass M core is traditionally used to discriminate CCSN from ECSN.Fryer et al. (2012) suggested the following criteria: if M core < 1.83 M ⊙ , there will be a degenerate CO core, which ends up as a CO WD; if M core > 2.25 M ⊙ , a non-degenerate CO core is formed, and subsequent thermonuclear burning leads to the formation of an FeNi core, which ultimately collapses to be a NS or BH; if 1.83 M ⊙ < M core < 2.25 M ⊙ , the star will form a partially degenerate CO core.If the CO core reaches a critical mass M CO,crit (= 1.08 M ⊙ ), it will non-explosively burn into a degenerate ONe core.If the ONe core mass could increase to M ecs = 1.38 M ⊙ , the star eventually collapses into a NS through an ECSN, or otherwise forms an ONe WD.
In addition, Case BB MT may strip most of the envelope of stars before SN explosions (Tauris et al. 2013(Tauris et al. , 2015(Tauris et al. , 2017)).For such ultra-stripped supernovae (USSN), we use the recipe in the the rapid model to estimate the compact star mass.
Taking into account all the above factors, we adopt the following five models for SN explosions: ( The natal kick velocities imparted on the NSs are set as follows. (1) In both rapid and delayed SN models, the kick velocity is assumed to follow a Maxwellian distribution with the dispersion velocity σ k = 150 or 300 km s −1 for CCSN, and σ k = 20 or 40 kms −1 for ECSN.

Deng et al.
(2) The remnant masses and natal kicks in stochastic SN satisfy some specific probability distributions depending on the masses of the CO core (Mandel & Müller 2020).
(3) For USSN, we use the formula given by Andrews & Zezas (2019) to evaluate the kick velocity (the best fit model of the SN kicks in Bray & Eldridge 2018): where M He and M NS are the masses of the pre-SN helium star and the NS, respectively.

Binary evolution models
Adopting different MT models, SN models, kick velocity prescriptions, and CE efficiency α, we construct 324 binary evolution models, and present the details in Table 2.The name of each model describes the properties of the key parameters.For example, model MTISNAk20150α0.1 means [MT Model I + (SN A) + (σ k = 20/150 kms −1 for EC/CCSN) + (α CE = 0.1)].We evolve 10 8 binary systems and obtain a few to several hundred thousand newborn DNS systems in each run.
Using the formulae of Peters (1964), we trace the orbital evolution of each DNS with time driven by GW radiation.We limit the evolutionary time to the Galactic age of 12 Gyrs, and calculate the total number of the Galactic DNSs.Each DNS spends specific time in the pulsar phase, and the number and characteristics of observable Pulsar-NS systems are dependent on the pulsar's evolution, the sensitivity of the radio telescopes, and selection effects.They are discussed in the following subsection.The first-born NS in DNSs behaves as radio pulsars mainly after the birth of the second NS.Because it had been (partially) recycled during the CE evolution and the subsequent Case BB MT phase, its spin period and magnetic field at the birth of the second NS (or at the birth of the PSR-NS binary) should deviate from its initial values.Three-dimensional hydrodynamic simulations by MacLeod & Ramirez-Ruiz (2015) show that NSs undergoing typical CE episodes could accrete a few percent their own mass or less depending on the density gradient inside the envelope and the assumption of NS's hypercritical accretion under neutrino-dominated cooling.The additional phase of MT from the helium star to the NS after the CE stage would also be responsible for the spin and magnetic field of the firstborn pulsar (Jiang et al. 2021, and references therein).However, the uncertainties in the accreted mass during both processes and the mechanism(s) of accretion-induced magnetic field decay make it difficult to evaluate their influence on the spin and magnetic field evolution in a very reliable way.We instead use the current spin periods and magnetic fields as references for the pulsars.The spin periods of the recycled pulsars in Galactic DNSs are distributed in the range of 17 − 185.5 ms (see Table 1).Assuming that the NSs are spinning down due to magnetic dipole radiation, from their spin periods and period derivatives we obtain the magnetic fields ∼ 10 8.5 − 10 10.3 G. Therefore, we tentatively assume that, at the birth of the second NS, the spin periods of the first-born NSs are uniformly distributed in 10 − 200 ms, and the initial magnetic fields (in the logarithm) are uniformly distributed in 8 − 11.We note these choices are somewhat arbitrary, and in section 4.3 we vary our assumptions of the spin periods and magnetic fields and discuss their possible influence.

Radio characteristics of DNSs
We also assume that, the first NS's surface magnetic field decays spontaneously in an exponential form  where B min is the lowest magnetic field strength taken to be 10 8 G (Os lowski et al. 2011), and τ d is the magnetic field decay timescale.We set it to be 10 9 yr according to the best-fitting models of Chattopadhyay et al. (2020) and Sgalletta et al. (2023).
We assume that the first-born NS acts as a pulsar and starts spinning down by magnetic dipole radiation after the second NS's birth.We adopt the fitting formula in Faucher- Giguère & Kaspi (2006) to estimate the pulsar's radio luminosity log L = log Here, L 0 = 0.18 mJy kpc −2 , Ṗ = (10 −15 ss −1 ) Ṗ15 is the period derivative, ϵ P = −1.5, ϵ Ṗ = 0.5, and L corr is a zero-centered normal distribution with σ Lcorr = 0.8.The flux of the pulsar is where D is the distance between the pulsar and the Solar system.The pulsar could be observed if F is larger than the limiting flux S min of the radio telescope.Approximately half of the Pulsar-NS systems in the Galaxy were detected by the Parkes Multi-beam Pulsar Survey.We set the limiting flux of the Parkes telescope to be 0.2 mJy (Manchester et al. 2001).
The distance between a pulsar and the Solar system depend on the initial location and the subsequent motion of the pulsar under the Galactic gravitational potential.Following Faucher-Giguère & Kaspi ( 2006) we set the birth distribution to be where r birth is the pulsar's birth distance from the Galactic center, R peak = 7.04 kpc and σ r birth = 1.83 kpc.Combining the Galactic rotation curve and the kick velocity at the second SN, one can follow the motion of the DNS.We use the NIGO code 3 to calculate the three-dimensional (3D) movement of the Pulsar-NS system and the Solar system with time, and to estimate their distance according to their 3D positions.The gravitational potential in the NIGO code includes contributions from the Galactic bulge, the disk, and the dark matter halo (Miyamoto & Nagai 1975;Flynn et al. 1996;Navarro et al. 1997;Prugniel & Simien 1997).The pulsar's radiation is known to be anisotropic.The beaming effect significantly influences the observable number of pulsars.We adopt the formula of Tauris & Manchester (1998) to model the beaming fraction, where P is the spin period of the pulsar in seconds.Finally, when the pulsars cross a "death line" in the P − Ṗ diagram, they stop emitting in the radio band (Rudak & Ritter 1994), log 10 Ṗ = 3.29 × log 10 P − 16.55, (13) log 10 Ṗ = 0.92 × log 10 P − 18.65.( 14) Here, we adopt the second equation ( 14).Following Sgalletta et al. (2023), we add a cut-off on the efficiency of radio emission as in Szary et al. (2014) to avoid the accumulation of pulsars at the death lines.The radio efficiency is defined as ξ = L/ Ė, where Ė is the pulsar's spin-down power.We set a threshold of ξ max = 0.01.If ξ > ξ max , the pulsar stops emitting radio waves (Szary et al. 2014).

RESULTS
3.1.The Orbital period and eccentricity distributions of Galactic DNSs

Qualitative analysis
3 Numerical Integrator of Galactic Orbits, which was designed to predict the orbital evolution of test particles moving within a fully-analytical gravitational potential generated by a multi-component galaxy.NIGO has been included in the Astrophysics Source Code Library and can be downloaded from https://github.com/LucaJRossi/NIGO(Rossi 2015).The gray symbols represent the DNS systems that need further confirmation.The black dash lines represent the merger time of tgwr = 10 8 , 10 9 , and 1.2 × 10 10 yrs, from left to right, respectively.Note that if the observable number is less than 0.01, the color of the pixel grid (∆P orb − ∆e) is set to white.
We plot the distribution of the Galactic DNSs in the P orb − e plane (Figs.1-4).Fig. 1 shows the impact of α CE on the number and orbit of DNSs in the SN A model.The number of observable DNSs increases by about 5 times when α CE increases from 0.1 to 10.The modeled DNSs are mostly distributed along the strip extending from lower bottom to upper right (i.e., from near circular short orbits to eccentric wide orbits).While quite a few observed DNSs are indeed on this strip, most of the model distributions face difficulty in accounting for DNSs with relatively long orbital periods and small eccentricities.
Fig. 2 shows the effect of kick velocities in the SN A model with α CE = 1.The larger kick velocity, the fewer DNSs formed (the observable DNS number in the model with k20150 is approximately four times that in the model with k40300).When the second SN is CCSN, most DNSs are short-period and high-eccentricity systems, and it is difficult (with σ k = 150 km s −1 ) or even impossible (with σ k = 300 km s −1 ) to form long-period systems.When the second SN is ECSN, most DNSs are low-eccentricity and long-period systems, but very few of them can merge within the Galactic age.In summary, both high and low kicks are required for the formation of the Galactic DNSs.

Quantitative analysis
We use the Bayesian analysis method employed by Andrews et al. (2015) to quantitatively compare the computed results with observations.We consider two independent orbital parameters P orb and e for each model.According to the Bayesian theory, the posterior probability related to model M given the observed data D is where D denotes P orb and e, P (D|M ) is the likelihood of observed data given that the model is true, and we set it to be Λ(D), P (M ) is the prior probability of a specific model, and P (D) is a normalizing constant which is independent of the model.Since all models have the same prior probability, P (M )/P (D) can be normalized as a constant C. Thus, P (M |D) is rewritten to be The relative probability of each model is determined by Λ(D) because C is independent of the models.Since the observed DNS systems are independent, the relative probability is equal to the probabilities of each independent system, Λ(D) = Λ (log P orb , e) Here P (log P orb,i , e i |M ) is the probability density of a model at a specific point in the log P orb − e parameter space.Similar as Chu et al. (2022), for each pixel (log P orb,i , e i ), we adopt a 2D normal distribution N (log P orb , e|log P orb,i , e i , σ 2 log P orb,i , σ 2 ei )| M to avoid zero probability at some pixels.The standard deviations of log P orb,i and e i are constant and assigned to be σ log P orb,i = 0.3 and σ ei = 0.03.So far, 18 DNSs have been observed in the Galactic field, of which 3 DNSs have not been fully confirmed (see Table 1).Therefore, we consider two sets of data.The first set of data D 1 consists of 15 fully confirmed DNSs, and the second set D 2 includes all 18 DNSs.Fig. 5 shows the values of Λ(D 1 ) as a function of the CE parameter α CE .Larger value of Λ(D) means that the model predictions are more consistent with observations.
We draw the following conclusions from Fig. 5: • The values of Λ(D) vary in a wide range with different input parameters.There is not apparent dependence of Λ(D) on the CE efficiency α CE and the kick velocity distribution.
• Overall nonconservative MT I and II models are considerably better than nearly conservative MT III models.
• SN B model (rapid CCSN explosion and M ecs = [1.83− 2.75] M ⊙ for ECSN) is more preferred compared with other SN models.

The merger rate of Galactic DNSs
Another way to calibrate the models of the formation and evolution of DNSs is comparing the expected DNS merger rate in the Galaxy by first principles with that inferred from the known Galactic population of Pulsar-NS systems.2015) derived the merger rate to be 21 +28 −14 Myr −1 at 95 per cent confidence.The method was recently applied by Pol et al. (2019Pol et al. ( , 2020) ) with updated DNS data, including the highly eccentric DNS system PSR J0509+3801 (Lynch et al. 2018), to estimate the merger rate to be 37 +24 −11 Myr −1 at 90 per cent confidence.Grunthal et al. (2021) revisited this problem taking account of recent observation of the longitudinal and latitudinal beam shape of PSR J1906+0746 (Desvignes et al. 2019), and derived a Galactic DNS merger rate of 32 +19 −9 Myr −1 at 90 per cent confidence.Fig. 6 compares the calculated and inferred DNS merger rates.For MT I and III models, to be compatible with these estimates requires α CE ∼ 5 − 10, while for MT II models, the allowed values of α CE reside in a wide range of 0.5 − 10.Among the ten best models selected by the Bayesian method in the Section 3.1.2,only models MTIISNBk20300α0.5,MTIISNBk40300α0.5, and MTIISNBk40300α10 match the derived DNS merger rate.So, we conclude that the best fit model is MTIISNBk20300α0.5.

INFLUENCE OF THE DOPPLER SHIFTING EFFECT
In the above section we have neglected the possible influence of the Doppler shifting effect and assumed that all Pulsar-NS binaries have the same detection probability.Actually it is challenging to detect eccentric binary pulsars with short orbital periods in pulsar searches considering the Doppler modulation of the pulsed signal as the pulsar moves with respect to the center of mass of the binary system during an observation (Johnston & Kulkarni 1991;Bagchi et al. 2013).Johnston & Kulkarni (1991) developed an analytical framework for calculating the reduction in signal-to-noise ratio caused by binary motion in the case of circular orbits.Bagchi et al. (2013) generalized the work of Johnston & Kulkarni (1991) and presented analytical expressions to calculate the sensitivity reduction for orbits of arbitrary eccentricity.Based on the results of Bagchi et al. (2013), Chattopadhyay et al. (2021) derived a fitting formula in the case of a 1.4 M ⊙ NS and an observation duration of 1000 s at an inclination angle of 60 • , and fit eccentricities via linear regression for e = 0.1, 0.5, and 0.8.These authors suggested that binary pulsars are detectable if where where m m = −8.9,c m = −27.68,m c = −3.4,and c c = −5.72.
The above criterion depends on both the orbital evolution and the pulsar's spin evolution.The latter is determined by the NS spin period and magnetic field when it behaves as a radio pulsar.In section 2.2 we simply adopt uniform distribution for both P s and log B inferred from their current observational data.However, due to previous accretion during the CE stage and the Case BB MT stage, both parameters should be related to the orbital period and eccentricity of the DNS systems4 .Moreover, the birth5 P orb and P s of DNS systems are not be exactly reflected by their currently observational data due to the loss of rotational energy and the decay of orbits caused by GW radiation.Taking these effects into account in a qualitative manner and placing particular weight on the very wide-orbit DNS PSR J1930−1852 (Swiggum et al. 2015), Tauris et al. (2017) obtained a simple empirical correlation between their birth P s and P orb , P s = (36 ± 14) ms P orb (day) 0.4±0.1 . (21) We adopt above equation to evaluate the birth spin periods of the first-born pulsars.
As to the pulsar's birth magnetic field, we select ten relatively young DNSs from Table 1 with the characteristic ages less than 10 9 yr, to ensure that their magnetic fields have not substantially decayed since the birth of the DNSs.Assuming that the characteristic age reflects the true age, we use the GW-induced orbital evolution functions and a Monte Carlo method to obtain the birth P orb and e.Then we employ a nonlinear least squares method to derive the correlation between the magnetic field and the orbital period and eccentricity (see Fig. 7), log B(G) = 1.47e 0.051 + 8.34P orb (day) 0.018 . ( The coefficient of determination of the above fit is 0.764.Note also that the results of Bagchi et al. (2013) suggest a stiff cut-off of radio pulsar detection, which means that no pulsar can be detected if Eq. ( 18) is not satisfied.We instead revise it to be a soft cut-off: if a binary crosses the cut-off line, the observable number is reduced by a factor According to the results of Bagchi et al. (2013), we find that f d is limited in the range of 0.3 − 0.55.Fig. 8 compares the cumulative distributions of the orbital period (P orb ), eccentricity (e), magnetic field strength (B), and spin period (P s ) (from top left to bottom right) of DNSs for both observational and modeled DNSs.We  find that considering a stiff cut-off due to the Doppler shifting effect significantly reduces the number of DNS systems with P orb < 1 day, inconsistent with observations, and the modeled distributions of e, P s and B also deviate from observations.On the other hand, the predicted orbital distribution of DNSs without considering the Doppler shifting effect is incompatible with observations.In comparison, when soft cut-off considered, the distributions of all the four parameters seem to better match observations.
We then re-calculate Λ(D 1 ) taking into account the Doppler shifting effect (see Table A), and present the results in Fig. 9. Comparing it with Fig. 5, we see the following features.
• Considering a soft cut-off improves most of the models in general.
• Similar as in Fig. 5, MT I and II models are better than MT III model, and SN B model is more preferred compared with other SN models.• In most cases, σ CCSN = 300 km s −1 and σ ECSN = 20 km s −1 are preferred over other kick velocity combinations.
In addition, we find that adopting the prescription proposed by Bray & Eldridge (2018) for the kick velocity of USSN (dark lines) does not significantly change the properties of the DNS population.However, assuming that USSN kicks follow the same distribution as ECSN kicks (i.e., σ k = 20/40 km s −1 , cyan lines) can significantly improve the match between SN E models and observation.
Combining with Fig. 6, we find that among the ten best models selected by the Bayesian method, model MTIISNBk20300α5 seems to best match the derived DNS merger rate.

DISCUSSION AND CONCLUSIONS
In this work, we construct 324 binary evolution models to investigate the formation of Galactic DNSs, considering the effects of MT stability, SN mechanisms and related kick velocity distribution, and CE efficiency.We evolve 10 8 binary systems and obtain a few to several hundred newborn DNS systems in each model.Using the formulae of Peters (1964), we trace the orbital evolution of each DNS with time driven by GW radiation.Combining with empirical relations for known Galactic pulsars, we predict the radio characteristics and merger rates of the DNSs in the Galaxy.We aim to evaluate the influence of the key parameters by comparing theory with observation.

CE evolution and MT stability
The evolution of the CE plays a key role in the formation of DNSs, but the specific physical processes are still not well understood (Ivanova et al. 2013).Regarding the formation of DNS in the Galaxy, different groups have used different CE parameters in their studies.Andrews et al. (2015) combined α CE with λ and set α CE λ = 0.01, 0.05, 0.1, 0.2, 0.25, 0.3, and 1.0.Using the Bayesian analysis method, they concluded that α CE λ = 0.3 − 1 is consistent with pulsar observations, while α CE λ ≤ 0.25 can be effectively ruled out, with the best value being 0.5.Considering that the value of λ for supergiants is less than 0.1 (Dewi & Tauris 2000;Podsiadlowski et al. 2003;Wang et al. 2016), their results indicate that α CE should be ≳ 5. Vigna-Gómez et al. (2018) used the values of λ calculated by Xu & Li (2010b,a) and set α CE = 0.1, 1 and 10 to simulate the formation of Galactic DNSs.Their results show that the value Deng et al.
of α CE mainly affects the birth rate and merger rate of DNSs, but the Bayesian analysis neither clearly favours nor disfavours the α CE variations with respect to the fiducial value 1. Kruckow et al. (2018) calculated the binding energy parameter λ of star's envelope based on their evolutionary state and set α CE = 0.2, 0.5, and 0.8, and showed that the birth rate and merger rate of DNSs increase with increasing α CE .Chu et al. (2022) took λ to be 0.5 and α CE = 0.1, 1, and 10.Their results of Bayesian analysis showed that the formation of Galactic DNSs requires a high efficiency factor (α CE = 10).Sgalletta et al. (2023) calculated the value of λ based on the model given by Claeys et al. (2014), and set α CE = 0.5, 1, 3, and 5. Comparing their results with the predicted DNS merger rate (37 +24 −11 Myr −1 in Pol et al. 2020), they concluded that only α CE = 3 matches the observation.
The occurrence of CE evolution is critically dependent on the stability of MT (Soberman et al. 1997;Shao & Li 2014;Mandel & Broekgaarden 2022).Packet (1981) pointed out that, when a star accretes about 10%∼15% of its original mass, it can be spun up to the critical rotation, and the centrifugal force would inhibit further accretion (Petrovic et al. 2005).However, in close binaries tidal synchronization can prevent the star reaching the critical rotation (de Mink et al. 2013).In addition, stars with strong magnetic fields may decelerate their rotation under the effect of magnetic braking (Dervişoǧlu et al. 2010;Deschamps et al. 2013;Shao & Li 2014).
Be/X-ray binaries are often used to evaluate the efficiency of stellar accretion (e.g., de Mink et al. 2013;Shao & Li 2014;Vinciguerra et al. 2020).Shao & Li (2014) investigated the binary interaction channel for the formation of Be/X-ray binaries.They showed that, if adopting the rotation-limited accretion model (Model I), the masses of the accreting stars are mostly less than 6 M ⊙ , which are inconsistent with the masses of the Be stars.Vinciguerra et al. (2020) compared the observed population of Be/X-ray binaries in the Small Magellanic Cloud with simulated populations of Be/X-ray binaries.Their results suggest that at least 30 per cent of the mass donated by the progenitor of the NS is typically accreted by the B-star companion.Additionally both Shao & Li (2014) and Vinciguerra et al. (2020) disfavor conservative MT.On the other hand, it has been suggested that the rotation-limited accretion model may help form the Galactic Wolf-Rayet/O binaries (Shao & Li 2016), Be/black hole binaries like MWC 656 with an orbital period of ∼ 60 days (Casares et al. 2014), and OB star-black hole system binaries (Langer et al. 2020), without the occurrence of CE evolution (Podsiadlowski et al. 2003;Shao & Li 2014;Grudzinska et al. 2015).
In this work, we adopt three MT models, the calculated λ values by Xu & Li (2010a,b) and Wang et al. (2016), and α CE = 0.1 − 10 to simulate the formation of DNS in the Galaxy.We find that the Pulsar-NS distribution does not effectively constrain the MT model and the value of α CE , but the derived DNS merger rate prefers MT II model and α CE = 0.5 ∼ 10.However, our best fit models with and without considering the Doppler shifting effect reveal α CE = 5 and 0.5, respectively.The former is broadly compatible with the studies on the formation of black hole low-mass X-ray binaries (e.g., Kalogera 1999;Podsiadlowski et al. 2003;Kiel & Hurley 2006;Yungelson & Lasota 2008), while the latter is consistent with studies on the formation of white dwarf binaries (e.g., Zorotovic et al. 2010;Davis et al. 2012;Nandez et al. 2015;Ge et al. 2022;Zorotovic & Schreiber 2022;Scherbak & Fuller 2022).Since the progenitors of black holes and white dwarfs are respectively mass stars and intermediate-/low-mass stars, the difference could imply possible dependence of α CE on stellar mass and mass ratio.

SN models and natal kicks
In this work, we adopt rapid, delay and stochastic CCSN models to estimate the CO masses.We note that under the same initial conditions, there are no significant difference between the predicted PSR-NS populations in the rapid and delay SN models, while the performance of the stochastic SN model seems not be as well as the other two.The main cause is associated with its different prescriptions of the kick velocity distribution.There is also large difference between the DNS merger rates calculated with the rapid/delay models and with the stochastic SN models.
Besides CCSN, we also consider the contribution from ECSN and USSN.We assume two possible helium core mass range M ecs = (1.83−2.25)M⊙ , and (1.83−2.75)M⊙ for ECSN.Our simulations show that using M ecs = (1.83−2.75)M⊙ better match the observations than using M ecs = (1.83− 2.25)M ⊙ , indicating that a large portion of NSs in DNSs have been formed via ECSN with a relatively small kick velocity (see also Shao & Li 2018a).
We adopt Maxwellian distributions for the kick velocities of CCSN and ECSN.For the specific combinations of the σ k values (20/150, 40/150, 20/300, 40/300 km s −1 ), our Bayesian analysis results indicate that σ k = 20/300 km s −1 is better than other combinations in most cases.Furthermore, if taking σ k = 300 km s −1 for CCSN, only MT II model can match the merger rate of Galactic DNSs.
As to USSN, it is not certain whether the NS is exerted a significantly smaller kick compared with CCSN.For example, Vigna-Gómez et al. (2018) adopted a Maxwellian distribution with σ k = 30 km s −1 for USSN kick velocities following Podsiadlowski et al. (2004), but numerical simulation by Schneider et al. (2021) suggested that stripped stars on average give rise to lower-mass NSs, higher explosion energies, and higher kick velocities (σ k = 315 ± 24 km s −1 ).In this work, we use the prescription proposed by Bray & Eldridge (2018) for the kick velocity of USSN, and find that there are not significantly changes in the properties of the DNS population and their merger rate.If we instead assume that USSN kicks follow the same distribution as ECSN kicks (i.e., σ k = 20/40 km s −1 ), then the match between SN E models and observation is greatly improved, as shown in Fig. 9 (in the cyan color and labeled as SN Es).Therefore, the USSN kick distribution will be an important subject for further research.

Summary
In this work we construct 324 models for the formation of Galactic DNSs taking into account various possible combinations of the key parameters.We try to explore their influence by comparing the properties of modeled DNS population with observation.Our work differs from previous studies in several ways.For example, compared with Andrews et al. (2015) and Chu et al. (2022), we used the λ values calculated from stellar evolution and considered the radio properties6 and merger rate of DNSs to constrain models.Our results show that the double-core CE channel (Channel II) also contribute to the formation of the Milky Way DNSs, which was not considered by Andrews et al. (2015).Compared with Vigna-Gómez et al. ( 2018) and Sgalletta et al. (2023), we considered more choices for the MT efficiency.We systematically explore the effects of different MT models, CE efficiency factor α CE , SNs and NS kicks on the formation of Galactic DNSs, and obtain useful estimates of the key parameters by comparing models predictions with both observed Galactic DNSs and their inferred merger rates.We list all the adopted models in Table 2 and use Bayesian analysis to evaluate them (see Section 3.2).Our main results are summarized as follows: 1. Comparing the predicted and derived Galactic DNSs merger rates reveals the allowed CE efficiency factor α CE ranging from 0.5 to 10.However, the Bayesian analysis indicates that there is no apparent dependence of the Pulsar-NS distribution on α CE .We find that α CE = 0.5 and 5 in our best fit models when the Doppler Shifting effect is not considered and when a soft cut-off related to the Doppler Shifting effect is considered, respectively.
2. Adopting MT II model (with 50% accretion efficiency) would allow more models to match the Galactic DNS merger rates, while Bayesian analysis indicates that there is no apparent dependence of the Pulsar-NS distribution on the MT mode.
3. Model SN B (with M ecs = [1.83− 2.75] M ⊙ for ECSN) seems to significantly better match observation than other SN models, implying that ECSN plays an important role in the formation of Galactic DNSs.
4. In most cases, σ ECSN = 20 km s −1 and σ CCSN = 300 km s −1 combination are preferred compared with other kick velocity combinations.The results are also sensitively dependent on the kick prescriptions for USSN.
Considering the fact that the current DNS sample is rather limited and there are many parameters involved in the formation of DNSs, we caution that our results are preliminary and subject to several issues.For example, we assume that the Milky Way has a constant metal abundance and a constant SFR, which will affect the merger rate of DNSs, even though the impact on the current merger rate is limited according to previous researches (Chu et al. 2022;Sgalletta et al. 2023).We do not take into account the merger history of the Milky Way when calculating the 3D movement of DNSs, which might have significant influence on the spatial distribution of DNSs.In addition, we assume that part of the accreting material is ejected out of the binary in the form of isotropic wind, while the actual situation for mass loss may be more complicated, which impacts the orbital evolution (Soberman et al. 1997;Kruckow et al. 2018;Vinciguerra et al. 2020).Together with more detections of DNSs in the Milky Way, a thorough investigation on the evolutionary history of DNSs needs to incorporate the formation of massive binaries (e.g., Wang et al. 2020), high-mass X-ray binaries (e.g., Misra et al. 2023), and binary pulsars in a self-consistent way.This will be the subject of the future study.

APPENDIX APPENDIX A: RESULTS OF BAYESIAN ANALYSIS
Table A. The posterior probabilities, ranking numbers, birth rates and merger rates of Galactic DNSs in all adopted models.The second and third sets correspond to the results without and with considering the Doppler shifting effect, respectively.
Faucher-Giguère & Kaspi (2006) performed Monte Carlo-based population synthesis of the birth properties of Galactic isolated pulsars, their time evolution, and their detection in the Parkes and Swinburne Multibeam surveys.They presented a population model that appears generally consistent with the observations, in which newborn NSs have normally distributed spin periods (with µ P = 300 ms and σ P = 150 ms), and log-normally distributed magnetic fields (with µ log(B/G) = 12.65 and σ log(B/G) = 0.55).More recently,Chattopadhyay et al. (2020) andSgalletta et al. (2023) considered several different initial parameter distributions to model the radio characteristics of Galactic DNSs.Their best fit models indicate that the NSs' initial spin periods and magnetic fields are uniformly distributed in the range of 10 − 100 ms and 10 10 − 10 13 G, respectively.

Figure 1 .
Figure 1.Distribution of the Galactic DNSs in the P orb − e plane.Each panel demonstrates a snapshot of the current DNS population in a specific model (MTISNAk40300α0.1,MTISNAk40300α0.5,MTISNAk40300α1, MTISNAk40300α2, MTISNAk40300α5, MTISNAk40300α10).Different colors indicate the predicted DNS numbers, as seen from the color bars on the right.The stars, triangles, and circles represent the Galactic field DNS systems with subpopulation (i), (ii), and (iii), respectively.The square, plus, and cross represent the GC DNSs PSRs J1807-2500, B2127+11C, and J0514-4002, respectively.The gray symbols represent the DNS systems that need further confirmation.The black dash lines represent the merger time of tgwr = 10 8 , 10 9 , and 1.2 × 10 10 yrs, from left to right, respectively.Note that if the observable number is less than 0.01, the color of the pixel grid (∆P orb − ∆e) is set to white.

Fig. 3
demonstrates the influence of the initial mass ratio on the formation of DNSs.Similar toShao & Li (2018b), both channels I and II contribute to the formation of DNSs.In the SN B model, ∼ 60% DNSs form through channel I (especially for near circular short orbit DNSs), while in the SN A model, the numbers of DNSs formed through channels I and II are roughly comparable.Fig. 4 compares the results in different SN and MT models.SN B model produces approximately two times DNSs than other models, and their distribution seems more consistent with observations, especially with MT II.The main reason is that in this model M ecs = (1.83− 2.75)M ⊙ , which substantially increases the occurrence rate of ECSN.The relatively low kicks during the SN reduces the disruption probability of the binaries, and helps form DNSs with low-eccentricity.SN A and C models show similar distributions.The DNSs formed in SN D model are mostly

Figure 2 .
Figure 2. Distributions of Galactic DNS binaries in the P orb − e plane under the SN A model with αCE = 1.The top to bottom panels correspond to models MTISNAk20150α1, MTISNAk20300α1, MTISNAk40150α1, and MTISNAk40300α1.The left, middle and right panels represent the distributions of total DNSs, DNSs with the second SN being CCSN and ECSN, respectively.

Figure 3 .
Figure 3. Distributions of Galactic DNS systems.The top and bottom panels correspond to models MTISNAk20150α1 and MTISNBk20150α1, respectively.The left, middle, and right panels represent the distributions of total DNSs, DNSs formed through channel I (high initial mass ratio), and channel II (initial mass ratio close to 1), respectively.

Figure 4 .
Figure 4. Distributions of Galactic DNS systems in different SN and MT models.The top to bottom rows correspond to SN A-E models, and the left, middle and right columns MT I-III models, respectively.

Figure 5 .
Figure5.The posterior probability for different models.The red, blue, yellow, green and black colors correspond to SNA-E models, respectively.The solid, dotted, dashed and dash-dotted lines represent the cases with σECSN,CCSN=20/150, 20/300, 40/150 and 40/300 km s −1 .Note that for SN D and E models, the kick velocity distribution is not a single Maxwell distribution (see Table2).The cyan color corresponds to SNEs.Based on SNE and assuming that the kicks of both USSN and ECSN follow the same distribution (i.e., σUSSN = 20/40 km s −1 ).The left, middle, and right panels represent MT I, II, and III models, respectively.The best fit model is marked by a red star.

Figure 6 .
Figure 6.Same as Figure 5 but for DNS merger rates.The orange and cyan horizontal lines and the shaded areas represent the merger rates of 37 +24 −11 Myr −1 (Pol et al. 2019, 2020) and 32 +19 −9 Myr −1 (Desvignes et al. 2019), respectively.Note that the two red hollow stars represent the best fit model obtained by adopting and not adopting the Doppler shifting effect.

Figure 6
Figure6depicts the DNS merger rates calculated for different models.They are in the range of ∼ 1 − 70 Myr −1 depending on the model parameters.In general, the merger rate increases with the CE efficiency α CE and decreases with the kick velocity.The reason is that larger α CE results in more DNSs, and larger kick velocity causes more binary systems to be disrupted.However, Figure6also shows that the DNS merger rate does not change monotonically with α CE .In some cases the merger rate decreases with increasing α CE .This is because larger α CE results in wider post-CE binaries, which are more easily disrupted after the SN.Additionally, the merger rate is critically dependent on the MT models.The overall merger rates in MT I and II models are higher than in MT III model.On the observational side,Kim et al. (2003) developed a Bayesian statistical method to derive the merger rate probability distribution by modelling the Galactic DNS population as well as the selection effects based on the observed properties of known binaries and survey characteristics.Based on this work,Kalogera et al. (2004) estimated the merger rate ∼ 90 Myr −1 , considering PSRs B1913+16, B1534+12, and J0737−3039A.Considering PSRs B1913+16, B1534+12, J0737−3039A, J1756−5521, and J1906+0746 with estimated beaming correction factors, O'Shaughnessy & Kim (2010) obtained the merger rate ∼ 60 Myr −1 .Taking into account the known DNS binaries with the best observational constraints, including both PSRs J0737−3039A and B, Kim et al. (2015) derived the merger rate to be 21 +28 −14 Myr −1 at 95 per cent confidence.The method was recently applied byPol et al. (2019Pol et al. ( , 2020) )  with updated DNS data, including the highly eccentric DNS system PSR J0509+3801(Lynch et al. 2018), to estimate the merger rate to be 37 +24 −11 Myr −1 at 90 per cent confidence.Grunthal et al. (2021) revisited this problem taking account of recent observation of the longitudinal and latitudinal beam shape of PSR J1906+0746(Desvignes et al. 2019), and derived a Galactic DNS merger rate of 32 +19 −9 Myr −1 at 90 per cent confidence.Fig.6compares the calculated and inferred DNS merger rates.For MT I and III models, to be compatible with these estimates requires α CE ∼ 5 − 10, while for MT II models, the allowed values of α CE reside in a wide range of 0.5 − 10.Among the ten best models selected by the Bayesian method in the Section 3.1.2,only models MTIISNBk20300α0.5,MTIISNBk40300α0.5, and MTIISNBk40300α10 match the derived DNS merger rate.So, we conclude that the best fit model is MTIISNBk20300α0.5.

Figure 7 .
Figure 7.The birth magnetic field (B0) as a function of birth orbital period (P orb,0 ) and birth eccentricity (e0) for selected first-born (recycled) NSs in DNS systems.The derived selected DNS birth data are plotted with blue points.The blue line represents a fit correlation (Eq.[22]).

Figure 8 .
Figure 8. Cumulative distributions of the DNS parameters P orb , e, B and Ps (from top left to bottom right) for the model of MTIISNBk20300α5.The blue line represents the observational data, The orange, green, and red lines represents the calculated DNS population that does not take into account the Doppler Shifting effect, adopts hard (Eq.[18]) and soft (Eq.[23]) cut-off, respectively.

Figure 9 .
Figure 9. Similar to Figure 5, but with the Doppler Shifting effect taken into account.

Table 1 .
Parameters of the detected Galactic DNSs.The spin period P , period derivative Ṗ , magnetic field B, orbital period P orb , eccentricity e, mass Mpsr of the pulsar, mass Mcom of the companion NS, the total mass M total , the characteristic age τc, and the merging time τgwr due to GW radiation.

Table 2 .
Parameter settings for binary stellar evolution models

Table A
continued on next pageTable A (continued)

Table A
continued on next pageTable A (continued) Table A continued on next page Table A (continued) Table A continued on next page