On the role of supernova kicks in the formation of Galactic double neutron star systems

In this work we focus on a group of Galactic double neutron star (DNS) systems with long orbital periods of $ \gtrsim 1$ day and low eccentricities of $\lesssim 0.4$. The feature of these orbital parameters is used to constrain the evolutionary processes of progenitor binaries and the supernova (SN) kicks of the second born NSs. Adopting that the mass transfer during primordial binary evolution is highly non-conservative (rotation-dependent), the formation of DNS systems involves a double helium star binary phase, the common envelope (CE) evolution initiates before the first NS formation. During the CE evolution the binary orbital energy is obviously larger when using a helium star rather than a NS to expel the donor envelope, this can help explain the formation of DNS systems with long periods. SN kicks at NS birth can lead to eccentric orbits and even the disruption of binary systems, the low eccentricities require that the DNSs receive a small natal kick at the second collapse. Compared with the overall distribution of orbital parameters for observed DNS binaries, we propose that the second born NSs in most DNS systems are subject to small natal kicks with the Maxwellian dispersion velocity of less than $ 80 \,\rm km\,s^{-1} $, which can provide some constraints on the SN explosion processes. The mass distribution of DNS binaries is also briefly discussed. We suggest that the rotation-dependent mass transfer mode and our results about SN kicks should be applied to massive binary evolution and population synthesis studies.


Introduction
More than 40 years has passed since the discovery of the first double neutron star (DNS) system PSR B1913+16 (Hulse & Taylor 1975), and there are now over a dozen known in our Galaxy (see Table 1). The orbital periods P orb of the DNS binaries are distributed in a wide range of ∼ 0.1 − 45 days, around half are close binaries (P orb < 0.5 day) and the other half wide binaries (P orb > 1 day). The orbital eccentricities e vary in a range of ∼ 0.1 − 0.8, and ∼ 80% of them have a relatively low value of e < 0.4. The components of the DNS systems generally have a similar mass of ∼ 1.2 − 1.4M except for PSR J0453+1559. The observed pulsars are usually mildly recycled, with the spin periods P spin in the range of 17 − 185 ms and the surface magnetic fields of the order 10 8 − 10 10 G.
Many previous works (e.g., Flannery & van den Heuvel 1975;Tutukov & Yungelson 1993;Lipunov et al. 1997;Dewi & Pols 2003;Ivanova et al. 2003;Tauris et al. 2017) have investigated the formation and evolution of DNSs which involve a variety of binary evolutionary stages. The canonical channel for the DNS formation has already been established (Bhattacharya & van den Heuvel 1991;Tauris & van den Heuvel 2006, for reviews). In a primordial binary with both components more massive than 8M , the primary star evolves to firstly fill its Roche lobe and transfer material to the secondary star. After that the primary becomes a helium star, which finally explodes with a supernova (SN) to be a NS. Then the secondary star evolves and transfers mass to the NS, and the binary becomes a high-mass X-ray binary (HMXB). Because of the large mass ratio, the mass transfer is dynamically unstable and binary system enters a common envelope (CE) phase (see recent review by Ivanova et al. 2013) . During the CE phase, the orbital energy of the NS is used to dispel the donor's envelope, and the orbital period is dramatically decreased. As the stripped donor (a helium star) finally becomes a NS through SN explosion, a DNS is formed if the binary is not disrupted. The subsequent evolution of the DNS system is driven by gravitational wave radiation, and its orbital decay can lead to a DNS merger event associated with a short gamma-ray burst (Abbott et al. 2017).
Modelling the DNS formation is subject to a number of uncertainties, especially those in the processes of CE evolution, SN explosions and related natal kicks. CE evolution plays a vital role in the orbital evolution (Webbink 1984). During SN explosions the natal kicks imparted to newborn NS can cause the orbit to be eccentric or even disrupt the binary. It has been shown that most of the massive binaries are disrupted if the first NS receives a sufficiently large kick, and very few systems can survive to become HMXBs (Belczynski & Ziolkowski 2009;Shao & Li 2014). From the orbital parameter distribution of the observed HMXBs, Pfahl et al. (2002) suggested that the NSs in wide and nearly circular binaries likely received small natal kicks, implying that they originated from electron-capture supernovae (ECS, Podsiadlowski et al. 2004;van den Heuvel 2004) rather than core-collapse supernovae (CCS). After the second SN explosion, the survived DNSs are generally in eccentric orbits. Based on the distributions of the DNS masses and their orbital parameters, Beniamini & Piran (2016) proposed that the second collapse in the majority of the observed DNS systems involves small mass ejection and a low kick velocity (see also Tauris et al. 2017).
Since the recent discoveries of gravitational wave transients by the advanced LIGO detector (e.g. Abbott et al. 2016Abbott et al. , 2017, a study of gravitational wave astronomy enters a golden age. Based on massive binary evolution, the rates of the merger events in the local universe can be theoretically estimated by a binary population synthesis method. However, the input parameters and physical assumptions should be carefully verified before the method is used to predict the merger rates in the local universe. Population synthesis studies (see e.g., Portegies Zwart & Yungelson 1998;Andrews et al. 2015;de Mink & Belczynski 2015;Belczynski et al. 2018;Chruslinska et al. 2018;Kruckow et al. 2018;Vigna-Gómez et al. 2018;Giacobbo & Mapelli 2018a,b;Mapelli & Giacobbo 2018) can be used to evaluate these uncertainties by adopting different input parameters to match the calculated distributions of the DNS population with those in our Galaxy. In this paper, we model the population of the Galactic DNSs. Compared to previous works, we have updated some important treatments of binary evolution, paying particular attention to the stability of mass transfer in primordial binaries, CE evolution, SN explosions and natal kicks. Our results may provide helpful constraints on the input physics of massive binary evolution. The remainder of this paper is organized as follows. In Section 2, we introduce the binary population synthesis method and the adopted assumptions. We present our calculated results and discussions in Section 3. We make a comparison with previous works in Section 4 and give a summary in Section 5.

Binary population synthesis method
To explore the distribution of the Galactic DNS systems, we use the population synthesis code BSE originally developed by Hurley et al. (2002). With BSE we can simulate the evolution of a large of binary stars with different initial parameters. The modelling of binary evolution involves detailed treatments of stellar evolution, tidal interactions, mass and angular momentum loss, and mass and angular momentum transfer. In the code some modifications have been made by Shao & Li (2014). Several important physical inputs are described as follows.

Mass transfer in the primordial binaries
Mass transfer begins in a primordial binary when the primary star evolves to fill its Roche lobe. The mass transfer efficiency (i.e. the fraction of the accreted matter by the secondary among the transferred matter) is a vital factor to determine the fate of the binary system. Rapid mass accretion can drive the secondary star to get out of thermal equilibrium and the responding expansion may also cause the secondary star to fill its own Roche lobe, leading to the formation of a contact binary (Nelson & Eggleton 2001). The contact binary will evolve to merge if the primary is a main sequence star, otherwise the binary system will enter a CE phase (see e.g., de Mink et al. 2013;Shao & Li 2014). Packet (1981) showed that only a small amount of accreted mass can accelerate a star to reach its critical rotation. When the secondary star is a rapid rotator, the mass transfer efficiency is subject to large uncertainty (Langer 1998;Petrovic et al. 2005;Stancliffe & Eldridge 2009). When modelling the evolution of the primordial binaries, Shao & Li (2014) built three mass transfer modes. In Mode I, the mass accretion rate onto the secondary star is assumed to be dependent on its rotating rate, equal to the mass transfer rate multiplied by a factor of 1 − Ω/Ω cr (where Ω and Ω cr are the angular velocity of the secondary and the corresponding critical value, respectively); in Mode II, the mass accretion efficiency is simply assumed to be 0.5; in Mode III, the mass accretion rate is assumed to be the mass transfer rate multiplied by a factor of min(10 τṀ τ KH2 , 1), where τṀ is the mass transfer timescale and τ KH2 the secondary s thermal timescale (Hurley et al. 2002). In the case of rapid mass accretion, the secondary star will get out of thermal equilibrium, expand and become overluminous. The thermal timescale of the secondary is greatly decreased, so the mass transfer process in Mode III is generally quasi-conservative, which is disfavored in massive binary evolution (Shao & Li 2014) and will not be considered in this work. In comparison the mass transfer efficiency in Mode I can be 0.2 for Case B binaries (Shao & Li 2016). The maximal initial mass ratios for avoiding the contact phase can reach as high as ∼ 6 in Mode I, and drops to ∼ 2 in Mode II and III (Shao & Li 2014).
In our calculation we adopt Mode I to treat the mass transfer in the primordial binaries, unless specified otherwise. The reasons are briefly described as follows. When modelling the evolution of massive binaries, Shao & Li (2016) demonstrated that the Galactic Wolf-Rayet/O binaries can be produced only if the mass transfer is highly non-conservative. It is also noted de Mink et al. (2007) discussed the mass transfer efficiency in massive close binaries and found a tendency that systems with long initial orbital periods evolve in a more non-conservative way than narrow systems, which might be caused by the coupling of spin-up and tidal interaction during the mass transfer phase.
Highly non-conservative mass transfer in Mode I can prevent the binary entering a CE phase and result in the formation of a long period system (see Figure 6 of Shao & Li 2014). MWC 656 is the first confirmed Be/black hole binary with an orbital period of ∼ 60 days (Casares et al. 2014). Assuming that the rapid rotating Be star has been accelerated by stable mass transfer during the primordial binary evolution, Shao & Li (2014) showed that non-conservative mass transfer is more likely to lead to the formation of MWC 656like binaries. Recently Kruckow et al. (2018) suggested that the binary AS 386 may be a candidate Be/black hole system, which contains a Be star of mass 6 − 8M in a ∼ 131 day orbit. Such a long period system is hardly produced if the progenitor binary has experienced CE evolution (Podsiadlowski et al. 2003;Grudzinska et al. 2015).
By comparing the fraction of runaway Be stars among the whole Be star population, Boubert & Evans (2018) concluded that Be stars could be explained by an origin of mass transfer in binaries, that is, most of the runway Be stars are the products of disrupted Be/NS binaries. For the identified runway Be stars (Berger & Gies 2001), the stellar spectral types can be late to be B8 (the corresponding stellar mass ∼ 3M ). This feature can be well produced by inefficient mass accretion which has spun up the secondary star without significant mass growth. Now that we assume that the mass transfer is highly inefficient, for the matter that is not accreted by the secondary star, we let it escape the binary system in the form of isotropic wind, taking away the specific angular momentum of the secondary star. In addition, we follow Hurley et al. (2002) to treat the process of rejuvenation, during which the secondary star appears even younger when fresh hydrogen is mixed into its growing convective core due to mass accretion.

Supernova explosion mechanisms and natal kicks
A NS can be formed through both ECS and CCS. In the BSE code, the helium core mass at the base of the asymptotic giant branch is used to set the limits for identifying various CO cores (Hurley et al. 2002). If the helium core mass is lower than 1.83M , the star will form a degenerate CO core and finally leave a CO white dwarf. Above this mass the star develops a partially degenerate CO core. If the core reaches a critical mass of 1.08M , it will non-explosively burn into a degenerate ONe core. Stars with an ONe core more massive than 1.38M are believed to collapse into NSs through ECS, otherwise form an ONe white dwarf. When the helium core is massive enough to form a non-degenerate CO core, stable nuclear burning will continue until a CCS occurs. However, investigations show that the trigger of ECS explosions is subject to big uncertainties, and mass transfer in binary systems makes it more complicated (see e.g., Nomoto 1984;Podsiadlowski et al. 2004;Woosley & Heger 2015;Poelarends et al. 2017). In the calculations, we consider several possible helium core mass ranges ∆M ecs (1.83 − 2.25M , 1.83 − 2.5M , 1.83 − 2.75M ), within which the star may eventually explode in an ECS explosion. The upper masses of these ranges can directly distinguish whether a star ends its life in an ECS or a CCS explosion. The larger the mass range, the more NSs formed through ECS.
We adopt the rapid SN model of Fryer et al. (2012) to deal with the mass of NSs formed from CCS, which varies in the range of ∼ 1.1 − 2M . The gravitational mass of ECS NSs is set to be 1.3M considering the fact that some baryonic mass of the 1.38M ONe core is converted into the binding energy of the NS (Timmes et al. 1996;Fryer et al. 2012). A newborn NS is assumed to be imparted by a natal kick, which may lead to the disruption of the binary system. Based on the observations of the pulsars' proper motions, Hobbs et al. (2005) demonstrated that the pulsar birth velocities can be described by a Maxwellian distribution with a dispersion velocity of σ = 265 km s −1 . More recently, Verbunt et al. (2017) argued that the birth velocities of pulsars can be better fitted by a bimodal distribution with σ 80 and 320 km s −1 , which may respectively be related to ECS and CCS. Analysing the nature of observed DNSs reveals that the second formed NSs in most systems are likely to receive low kicks (Beniamini & Piran 2016;Tauris et al. 2017). In our calculations, the NS kick velocities are assumed to obey Maxwellian distributions. We adopt σ k,ecs = 20, 40 and 80 km s −1 for ECS NSs and σ k,ccs = 150 and 300 km s −1 for CCS NSs to examine their possible influence on the DNS formation.

Mass transfer in NS binaries
After the birth of the first NS, the binary evolves into the HMXB phase. Modelling the evolution of NS HMXBs reveals that the maximal mass ratio for stable mass transfer is ∼ 3.5 Podsiadlowski et al. 2002;Shao & Li 2012). If the donor is a helium star, we assume that the binary enters CE evolution only when the orbital period is less than 0.06 day (Tauris et al. 2015). In the case of unstable mass transfer that leads to CE evolution, we calculate the orbital energy of the binary system to examine whether it is sufficient to expel the stellar envelope. We use the results of Xu & Li (2010) and Wang et al. (2016) for the binding energy parameter λ of the envelope, and take the CE efficiency α CE to be 1.0 in the calculations.
For stellar wind mass loss, we employ the fitting formula of Nieuwenhuijzen & de Jager (1990), except for hot OB stars, for which we use the simulated rates of Vink et al. (2001).
For the stripped helium stars and Wolf-Rayet stars, we adopt the rates of Vink (2017) 1 . During the CE phase, we ignore any mass accretion because of its very short duration. In NS binaries with a helium donor star, mass transfer proceeds rapidly for a short (< 10 5 yr) time through Roche lobe overflow (Tauris et al. 2015), the accretion rate of the NS is assumed to be limited by its Eddington rate.

Results and discussions
We have performed a series of Monte Carlo simulations for a large number of binary systems with different initial parameters. For each of the considered models, we simulate the evolution of 10 8 binary systems. We follow Hurley et al. (2002) to set the initial parameters of the primordial binaries. The primary masses obey the initial mass function given by Kroupa et al. (1993), and the mass ratios of the secondary to the primary satisfy a flat distribution between 0 and 1. The distribution of the orbital separations is assumed to be flat in the logarithm in the range of 3 − 10 4 R . All binaries are assumed to have circular orbits, since the outcome of the interactions of binary systems with the same semilatus rectum is almost independent of eccentricity (Hurley et al. 2002). Several authors (Smith et al. 1978;Diehl et al. 2006;Robitaille & Whitney 2010) have proposed that the star formation rate in our Galaxy is ∼ 1 − 5M yr −1 . Here we adopt a moderate rate of 3M yr −1 over a period of 10 Gyr. When investigating the star formation history of the Galaxy, Rocha-Pinto et al. (2000) showed that the star formation rate is not likely to change by much more than a factor of two over the past 10 Gyr. So adopting a constant star formation rate should be an acceptable approximation.
From each of the population synthesis calculations, we obtain a few to several ten thousand initial DNS systems. The corresponding birthrates can be calculated according to the parameters of their primordial binaries (Hurley et al. 2002). After the DNS formation, we adopt the formulae of Peters (1964) to trace the subsequent orbital evolution with time due to gravitational wave radiation. The DNS merger rate within the lifetime of the Galaxy is mainly dependent on the initial orbital parameter of the DNS binary 2 . In the calculations, we assume that the DNS binaries evolve to merge when the orbital periods shrink to be less than 0.001 day. Since we adopt a constant star formation rate of 3M yr −1 , the birthrate of the Galactic DNS systems does not vary with time, while the merger rate tends to gradually increase over the past 10 Gyr. The total number of the DNS binaries in our Galaxy can be obtained by accumulating the product of the birthrate and the evolutionary time for every binary. In Table 2, we show the merger rate R merger , the birthrate R birth and the total number of the Galactic DNSs at present, and the fraction f of DNS systems with the second NS being formed through ECS among all DNSs. Obviously the larger values of both σ k,ecs and σ k,ccs result in the decreases of R merger and R birth . We obtain R merger ∼ 2 − 10 Myr −1 , R birth ∼ 3 − 27 Myr −1 , and the total number of DNSs in our Galaxy is of the order 10 4 . It is also seen that wider ∆M ecs leads to higher f , which varies in the range of 0.06 − 0.85.

The P orb − e diagram
Since the orbital evolution of DNS systems is solely driven by gravitational wave radiation, we can synthesize the orbital parameters of the simulated DNS population into a diagram that represents the snapshot after 10 Gyr star formation and evolution. We show in Figure 1 the calculated distribution of Galactic DNS systems (with the numbers given in Table 2) in the P orb − e plane when adopting ∆M ecs = (1.83 − 2.25M ). The left, middle and right panels correspond to σ k,ecs = 20, 40 and 80 km s −1 for ECS NSs, and the top and bottom panels correspond to σ k,ccs = 150 and 300 km s −1 for CCS NSs, respectively. Each panel includes 20 × 20 image elements, in which P orb varies logarithmically from 0.01 to 10 4 days with steps of 0.3, and e linearly from 0 to 1 with steps of 0.05. Until now there exist 14 observed DNSs whose orbital periods and eccentricities are precisely measured (see Table  1). They are plotted as blue triangles in the figure. Figure 1 shows that the calculated DNSs are mainly distributed in a region with relatively large eccentricities (e > 0.4), so the observed systems with long orbital periods ( 1 day) and low eccentricities ( 0.4) can be hardly reproduced. The reason is that with ∆M ecs = (1.83 − 2.25M ) the second NSs in most systems are born in a CCS explosion with a high natal kick, which tends to increase the orbital eccentricities and even disrupt the original systems (especially for wide binaries). Figures 2 and 3 show the similar distributions of Galactic DNSs in the P orb − e plane with ∆M ecs = (1.83 − 2.5M ) and (1.83 − 2.75M ), respectively. In these two cases the value of f is remarkably increased and can reach as high as 0.85 when σ k,ecs = 20 km s −1 and σ k,ccs = 300 km s −1 (see Table 2). The second born NSs in most systems are imparted to small natal kicks, and the long-period binaries can effectively avoid disruption. In addition, the predicted DNS binaries tend to have low eccentricities. Generally the observed DNS binaries can be well covered by a high probability region in the P orb − e diagrams when taking ∆M ecs = (1.83 − 2.75M ), expect the case with σ k,ecs = 80 km s −1 (see right panels). Consistent with previous works (e.g. Beniamini & Piran 2016;Tauris et al. 2017), our results suggest that the formation of the second NS be accompanied with a low natal kick.
We then try to quantitatively compare the calculated results with observations to constrain the adopted models with different input parameters. Although the observed Galactic DNS sample is limited by small size and the detection of such systems may be subject to several selection effects (see Section 2.1 of Tauris et al. 2017), we employ the Bayesian analysis used by Andrews et al. (2015) to compare the goodness of match for every model, which involves two independent orbital parameters of P orb and e. The Bayes' theorem gives a probability expression of where P (M | D) is the posterior probability, i.e. the probability related to our model M given the observed data D (i.e., the data set of P orb and e), P (D | M ) is the likelihood of the observed data given that the model is true, which is denoted as Λ(D), P (M ) is the prior probability of our model, and P (D) a normalizing constant that is independent of the model. As all models have the same prior probability, it can be absorbed into the normalizing constant as C = P (M )/P (D). So P (M | D) can be rewritten as Since C is independent of the models, the value of Λ(D) gives the relative probability of each model. Because all of the observed DNSs in the sample are independent, the probability of the observed data given the model is equal to the product of the probabilities of each binary system: Λ(P orb , e) = i P (P orb,i , e i | M ).
Note that the errors of the observed orbital period and eccentricity for each DNS system are very small and can be neglected. Then P (P orb,i , e i | M ) reflects the probability density of a model at a specific point in the P orb − e parameter space, which corresponds to a specific DNS. Adopting similar treatment of Andrews et al. (2015), we calculate P (P orb,i , e i | M ) by binning the simulated systems into two dimensions with a bin size of ∆ log(P orb ) = 0.3 in orbital period and ∆e = 0.05 in eccentricity. Finally the derived value of Λ(P orb , e) is used to rank our models (see Table 2). The model that best matches the observations involves the input parameters of σ k,ecs = 40 km s −1 , σ k,ccs = 300 km s −1 and ∆M ecs = (1.83 − 2.75M ).
To display the effect of the mass transfer process during the primordial binary evolution on the DNS population, in Fig. 4 we present the corresponding distributions of DNSs in the P orb − e plane when the mass transfer efficiency is taken to be a constant value of 0.5 (i.e., Mode II). The left, middle and right panels correspond to ∆M ecs = (1.83 − 2.25M ), (1.83 − 2.5M ) and (1.83 − 2.75M ), respectively. The kick velocity parameters are all set to be σ k,ecs = 40 km s −1 and σ k,ccs = 300 km s −1 . Compared with the situations in Mode I, the predicted distribution of the DNS systems shifts to a region with shorter orbital periods and seems not to match observations. The main reason is that the birthrate of close DNSs in Mode II is greatly enhanced since a large number of original NS/white dwarf binaries (post-CE products of HMXBs) in Mode I are transformed into DNS systems due to the higher mass transfer efficiency. Further, more efficient mass transfer can significantly increase the secondary mass, so it takes more kinetic energy of the NS to expel the secondary's envelope during the CE evolution.

The mass distribution of DNS systems
All of the observed Galactic DNSs except PSR J0453+1559 have similar masses of ∼ 1.2 − 1.4M (see Table 1). As expected, recycled pulsars are generally more massive than their non-recycled companions because a certain amount of mass is required to recycle a pulsar. In Figure 5 we present the distributions of the pulsar mass M p (first formed NS, left), the companion mass M c (second formed NS, middle) and the total mass M t (right) for the Galactic DNS systems under different assumed models. The panels form top to bottom correspond to ∆M ecs = (1.83 − 2.25M ), (1.83 − 2.5M ) and (1.83 − 2.75M ), respectively. In each panel, the black, red, and blue curves correspond to σ k,ecs = 20, 40, and 80 km s −1 , and the dashed and solid curves correspond to σ k,ccs = 150 and 300 km s −1 respectively. The green curves represent the observed data multiplied by a factor of 10 4 . We find that the calculated mass distributions of both DNS components generally have two peaks of ∼ 1.1M (minimum NS mass in the rapid SN model of Fryer et al. 2012) and ∼ 1.3M (the mass of NSs formed from ECS), in agreement with the results of Giacobbo & Mapelli (2018b) and Mapelli & Giacobbo (2018).
The calculated mass distribution at the peaks of ∼ 1.1M and ∼ 1.3M could be used to distinguish the NSs formed from CCS and ECS, respectively. It is unsurprising that the percentage at each peak changes with ∆M ecs , the former decreases and the latter increases with increasing ∆M ecs since more ECS NSs are produced. Even in the case of ∆M ecs = (1.83 − 2.25M ), quite a fraction of DNS binaries undergo at least one ECS explosion. This is consistent with the result of Giacobbo & Mapelli (2018a) for merging DNSs, indicating that the ECS is a fundamental process for the DNS formation. We can see that the overall distributions of calculated DNS masses, although in agreement with other studies, are actually not compatible with observations. The reason may be that the adopted SN mechanisms are still too simplistic, and there may not be one to one relation between the compact star mass and the core mass of their progenitors. Furthermore a NS can increase its mass due to the possible accretion during the recycling phase.
In spite of the mass distribution for DNS binaries, comparison between the calculated and measured orbital parameters gives the best choice of the input parameters of σ k,ecs = 40 km s −1 , σ k,ccs = 300 km s −1 and ∆M ecs = (1.83 − 2.75M ) (see the rank of different models in Table 2), which are therefore set to be default input parameters in the following.

Formation channels of DNS systems
The formation channels of DNS systems have been explored by many authors. Besides the canonical channel in which the binary experiences a HMXB stage (Bhattacharya & van den Heuvel 1991;Tauris & van den Heuvel 2006), there is an alternative channel which involves a double helium star phase. In the double core CE scenario (Brown 1995) a DNS system is the descendant of a double helium star binary, which is the result of a CE phase in a binary with both components of very similar masses. Both components of the binary have developed a compact core at the beginning of mass transfer. Such a binary will enter the contact phase that is followed by the CE evolution. The post-CE product is a double helium star binary in a relative close orbit. After experiencing two SN explosions, the survived binary is a DNS system (see also Dewi et al. 2006).
By analysing the calculated results, we find that the double helium star binary can be produced in an another way. Here both of the binary components also initially have almost equal mass. The secondary star is still in the main sequence at the onset of mass transfer. The mass transfer proceeds stably between the binary components, after which the primary leaves a helium star and the secondary approaches the end of main sequence. The secondary then evolves to be a supergiant star and transfers mass to the helium star, resulting in CE evolution. The outcome of CE evolution is a double helium star binary. A similar scenario was suggested by Portegies Zwart & Yungelson (1998), in which the mass exchange during the primordial binary evolution was assumed to be almost conservative, but the process of rejuvenation was not considered in their calculations. In the case of conservative mass transfer, this can cause the secondary to be younger, and stay in the main sequence when the primary explodes to be a NS. A highly non-conservative mass transfer like in Mode I can effectively reduce the rejuvenation of the secondary. The primary helium star firstly explodes to be a NS, which is subsequently recycled by the mass transfer from its helium star companion. When the secondary evolves to be a NS, the binary finally turns into a DNS system. Adopting the default model with the input parameters of σ k,ecs = 40 km s −1 , σ k,ccs = 300 km s −1 and ∆M ecs = (1.83 − 2.75M ), we obtain the total birthrate of Galactic DNSs to be ∼ 10 Myr −1 (see Table 2). The DNS binaries formed from the canonical channel have a birthrate of ∼ 4 Myr −1 , and the two scenarios involving a double helium star binary contribute a similar birthrate of ∼ 3 Myr −1 . This means that the double helium star channel is likely to dominate the formation of Galactic DNS systems. Figure 6 shows the calculated distribution of DNS systems in the P orb − e plane when taking into account different formation channels. The left and right panels correspond to the results from the canonical channel and the double helium star channel, respectively. In the canonical channel, the DNS systems are dominated by binaries with short orbital periods of 1 day and low eccentricities of 0.4, their distributions show similar features to that in Mode II (see the right panel of Figure 4). The double helium star channel tends to produce the DNS systems with long periods of 1 day and high eccentricities 0.1. The difference originates from the fact that the binary orbital energy that is used to expel the CE in the former is considerably smaller than in the latter. Here the fraction of systems that formed from the canonical channel among all DNSs is ∼ 0.3.
In Figure 7 we present the calculated distributions of the DNS systems in the P orb − e plane, in which the second born NS originates from ECS (left) or CCS (right) explosion. The NSs formed through CCS have received relatively large natal kicks, which tend to produce the DNS systems with large eccentricities. On the contrary, the DNS systems generally have low eccentricities for ECS explosions. A large number of relatively wide systems can survive the SN explosions, and evolve to DNS systems with long orbital periods and low eccentricities. The fraction of the DNS systems that experiences a second ECS explosion can be as high as 0.74 (see Table 2). We note that the PSRs J1757-1854 and B1913+16 have relative short periods and large eccentricities, and their orbital decay timescales are of the order 100 Myr. These two binaries are not easily accounted for by the current scenarios.

Comparison with other works
In this section, we compare our results with recent population synthesis works of Andrews et al. (2015), Kruckow et al. (2018) and Vigna-Gómez et al. (2018).
The most important parameters in Andrews et al. (2015) include the combined CE parameter α CE λ, the NS kick velocity distribution and the ECS mass range (same as our ∆M ecs ). A significant difference between their work and ours is the treatment of the α CE λ, which was set to be constant by Andrews et al. (2015). Note also that the double helium star channel is not considered by Andrews et al. (2015), who focus on the canonical channel. To match the orbital parameters of the known eight DNSs, their models are required to produce wider binaries than our canonical channel. This indicates that their adopted α CE λ is generally larger than ours. By testing different models to constrain the input parameters, they found that α CE λ 0.25 is effectively ruled out, the allowed range is 0.3 − 1.0 and the best choice is 0.5. Considering the fact that the value of λ for massive supergiants are less than 0.1 (Dewi & Tauris 2000;Podsiadlowski et al. 2003;Wang et al. 2016), this implies that the value of α CE in Andrews et al. (2015) is significantly larger than unity. After the discovery of PSR J1930-1852 with the longest known orbital period of 45 days, Swiggum et al. (2015) argued that many models of Andrews et al. (2015) cannot produce such wide DNS systems, probably requiring reconsideration of the binary evolutionary processes.
The investigation of Kruckow et al. (2018) involves the amplitudes of the NS natal kicks and the treatment of mass and angular momentum transfer. They proposed that a low kick velocity is applied when the envelope of the progenitor star is stripped. This can help produce the observed DNS systems with low eccentricities. Following the method of Soberman et al. (1997) to model the binary orbital evolution with Roche lobe overflow, Kruckow et al. (2018) assumed that the fractions of the mass escaped from the donor and ejected from the accretor are respectively 0.2 and 0.75, and the rest material (with the fraction of 0.05) is accreted by the accretor. In their default model the mass accretion efficiency is also very low, but the channel involving a double helium star binary was not included. For their modelled DNS systems, the probability distribution of the orbital periods is peaked at ∼ 1 − 3 days, slightly longer than that ( 1 day) predicted in our canonical channel. The reason is perhaps that a large fraction of the lost matter from the donor through a fast wind tend to increase the orbital separations of the progenitor systems (this can lead to more expansion of the donor and decrease the binding energy of the envelope during CE evolution, broadening the subsequent DNS systems). From their orbital parameter distribution (Figure 11), there is a shortage of wide systems with orbital period larger than several days with respect to the observed DNSs. Meanwhile close systems with orbital period less than ∼ 1 day are also deficient. Vigna-Gómez et al. (2018) suggested that a bimodal natal kick distribution (with σ k = 30 km s −1 and 265 km s −1 ) is preferred over a high-kick component alone. This is in accord with our assumption about the NS kick velocity distribution. They adopted the same treatments of the NS mass and kick velocities for both ultra-stripped SNe (see Tauris et al. 2015) and ECSs, which can boost the formation of DNS binaries with low eccentricities. In their fiducial model, the canonical channel is responsible for formation of about 70% of all DNSs and the channel with a double core CE phase approximately contributes about 21%. It can be seen that their double helium star channel contributes a significant fraction of DNS population, but still much smaller than that in our results. Since the probability distribution of the orbital parameters for DNS systems is not clearly shown by Vigna-Gómez et al. (2018), we cannot make further comparisons.

Summary
We have performed a population synthesis study of the formation of DNS systems in the Galaxy. Based on previous investigations on Wolf-Rayet/O binaries, Be/black hole systems and runaway Be stars, we argue that mass transfer in the primordial binaries is likely to be highly inefficient. So we adopt the rotation-dependent mass transfer mode (Mode I) in our calculations. We then investigate the formation channels of the DNSs and the effect of the SN kicks at the NS formation on the distribution of DNS systems. By comparing with observed DNS binaries, we can set important constraints on the model parameters. We summarise our main results as follows.
1. The second born NSs in most DNS systems have low kick velocities, with the dispersion velocity of less than 80 km s −1 in a Maxwellian distribution. This also requires that the helium core mass range for ECS is 1.83 − 2.75M if low kicks for ultra-stripped SNe are not taken into account.
2. The DNS population can be divided into two groups by the orbital period of ∼ 1 day, the systems with short periods mainly form through the canonical channel involving a HMXB phase and the wide binaries experience a double helium star binary stage during the progenitor evolution.
3. Our best model shows that the birthrate and merger rate of Galactic DNS systems in the canonical channel are ∼ 4 Myr −1 and ∼ 2 Myr −1 , and the corresponding rates in the double helium star channel are ∼ 6 Myr −1 and ∼ 2 Myr −1 , respectively.
At last we briefly discuss the applicability of the ∆M ecs that is ascribed to the process of ECS explosions. Although Andrews et al. (2015) extended the upper mass limit to be 3.24M in their study on DNS formation, our best range of ∆M ecs = (1.83 − 2.75M ) seems to be still wider than previously expected. Investigations indicate that, besides ECS, both ultra-stripped SNe and low-mass Fe-core SNe are expected to have considerably lower NS kicks than CCS of massive stellar cores (Tauris et al. 2015(Tauris et al. , 2017Janka 2017;Gessner & Janka 2018). This implies that the wide core mass range may include the contributions of ultra-stripped SNe and low-mass Fe-core SNe.
We thank the referee for constructive comments that helped improve this paper. This    In each panel, the black, red and blue curves correspond to σ k,ecs = 20, 40 and 80 km s −1 , and the dashed and solid curves correspond to σ k,ccs = 150 and 300 km s −1 respectively. The green curves denote the observed data multiplied by a factor of 10 4 . There are 7 (13) DNS binaries with component (total) masses precisely measured (see Table 1).