EVOLUTIONARY CHANNELS FOR THE FORMATION OF DOUBLE NEUTRON STARS

To simulate DNS formation, we adopt the population synthesis code StarTrack (Belczynski et al. 2008). We use a reference model with a given set of parameters, and we then vary a large number of them from the reference assumptions to assess their effect on DNS properties and formation channels.

StarTrack adapts the fitting formulae by Hurley et al. (2000) for single star modeling (time evolution of key macroscopic physical properties) and it modifies them to account for the effects of mass transfer to allow for modeling of binary evolution. The mass-transfer treatment (outside of the case of CE phases) has been tested against and calibrated to mass-transfer sequences (see Belczynski et al. 2008, and references therein). This code has been used for a wide range of studies and applications involving modeling of interacting binaries with all types of compact objects, and has been compared to and demonstrated agreement with observations of such binaries both in our Galaxy and in external galaxies since the early 2000's when the first version of the code was produced. It has also been maintained and kept up to date as our understanding of compact object formation and stellar winds has evolved over the years. In particular, updated prescriptions for massive stellar winds and their dependence on metallicity have been included (Belczynski et al. 2010). Although a recent revision of the code uses a compact object remnant mass function physically motivated by calculations of the SN mechanism (Fryer et al. 2012), our results rely on the remnant mass function calculated by Timmes et al. (1996). An additional recent parametrization of the envelope binding energy parameter λ, a measure of the central concentration of a stellar envelope which is relevant for CE calculations, was added to Startrack after our simulations were completed (Dominik et al. 2012) and is not included in the version of the code used here.

Initial conditions for each population model are as follows: the primary star in the binary is formed on the main sequence with a mass drawn from a power law distribution ∼M^−2.7 (Scalo 1986), while the secondary is determined based on the mass ratio value drawn randomly from a flat distribution between 0 and 1; the initial orbital separation is logarithmically flat and the initial eccentricity follows the dynamical-equilibrium distribution ∼2e. Allowing for wind mass loss, the primary star is then followed through the Hertzsprung gap (HG), red giant branch, He-main sequence, and asymptotic giant branch (AGB). If the primary is massive enough, it will collapse through a SN, leaving a compact remnant. Assuming the binary is not disrupted during the first SN, the secondary evolves through its lifetime, with possible interruptions by mass transfer sequences and also driven by wind mass loss and/or tides. Once it has evolved completely, if massive enough, the secondary will also leave a compact remnant as the product of a SN. Asymmetries during the SN impart an (assumed instantaneous with respect to the orbital period) kick to the NS or black hole, altering the shape of the binary orbit, following the equations of Kalogera (1996, 2000). The formation of a DNS requires that both stars have pre-SN masses within a specific range, such that the binary neither disrupts due to the SN kick nor does it merge prior to the second SN. The particular variables and prescriptions used are given in Section 2.1.1. The DNS population is given a constant birth rate over the past 10 Gyr. After the second SN, the orbital evolution due to gravitational radiation is followed using the equations of Peters (1964) until the current epoch.

2.1.1. Reference Model

2.1.2. Varied Parameters

In the standard DNS formation model, the most important phases of evolution occur after the primary becomes an NS. Once the secondary evolves beyond the main sequence, it may enter a phase of dynamically unstable mass transfer. As mentioned in Section 1, we completely ignore non-interacting systems because the short orbital periods in observed DNSs require mass transfer. The secondary will evolve onto the helium main sequence, possibly filling its Roche lobe upon evolving off the helium main sequence in the so-called Case BB mass transfer (Delgado & Thomas 1981). After forming a degenerate core, the star will collapse forming an NS in either an Fe core-collapse SN or an ECS.

Variations from the standard paradigm have been discussed in the literature or seen in population synthesis studies (Bhattacharya & van den Heuvel 1991; Brown 1995; Portegies Zwart & Yungelson 1998; Belczynski et al. 2002). We find the scenario outlined by Brown (1995) in which the initial mass ratio is close to unity, and the individual stars in the binary evolve off the main sequence together, forming double CEs occurs in less than 1 percent of systems, and ignore them in our analysis, (however see Schwab et al. 2010). While we find some NSs form through the accretion induced collapse of a WD the population is insignificant and we disregard them. Taking the standard formation paradigm outlined above, we can break down our models into three different evolutionary tracks depending on whether or not the system went through stable mass transfer after a CE (Case BB mass transfer), and, if it did, whether the second NS was formed in an ECS or an Fe core-collapse SN.

Channel I. As the simplest of our three formation channels, the systems that are created this way have the fewest constraints on their evolution. Therefore this DNS population covers the widest range in orbital period and eccentricity, which can be seen in the second panel in Figure 1. Following the formation of the first NS, the secondary evolves off the main sequence; the orbits of this population are wide enough that the massive star fills its Roche lobe when a significant convective envelope has developed and a CE ensues. The NS companion is massive enough that it eventually forms an Fe core before collapsing into the second NS in the system.

Channel II. In this channel, the binary orbit is tight enough following the CE event that once core He burning has finished the star expands, filling its Roche lobe again in Case BB mass transfer; provided the He star has not developed a fully convective envelope, mass transfer onto the first NS is stable (see Ivanova et al. 2003). The secondary evolves to form an Fe core and eventually collapses in a SN. Comparing panel two to three in Figure 1 shows that this stable mass transfer phase produces a DNS population with shorter orbital periods than that of Channel I.

Channel III. These systems follow a similar evolution to Channel II, except prior to core-collapse, the secondaries have lower He core masses, leading to NS formation through an ECS. These lower mass He star progenitors always expand significantly leading to the second phase of mass transfer onto the primary NS. Unless the system has a very small orbital separation, this phase of mass transfer is stable (Dewi et al. 2002; Ivanova et al. 2003). Since ECS kick velocities are small, the resulting DNS orbital geometries are more dependent on the mass loss of the secondary during NS formation. Due to the tight constraints on the mass of the secondary immediately prior to the SN (the secondary star has a mass within a narrow range such that carbon is allowed to burn, but oxygen is not), the resulting population has a very narrow distribution of geometries. This can be seen by the relatively narrow distribution of systems in Figure 1. Combined with the small kick velocities, systems going through this evolutionary channel tend to have systematically lower eccentricities.

A few known DNS systems have such characteristics that a number of different groups have been able to derive firm constraints on their evolutionary history, and we are able to confidently place them in the channels described here. Having this knowledge allows us to take a step further and, for the first time, impose evolutionary constraints on the population formation. In what follows we summarize the origin of these evolutionary constraints on the three best studied DNS systems. Once these constraints are imposed, we use the best models to draw conclusions about the most likely channels through which the rest of the DNS systems have evolved (Section 4.3).

J0737−3039. The binary pulsar is probably the most widely studied DNS in terms of its evolutionary history, which is now generally agreed upon. After the primary star became an NS, the secondary had to evolve off its main sequence filling its Roche lobe in the process and forming a CE. Both Dewi & van den Heuvel (2004) and Willems et al. (2004) conclude, based on the current geometry of the system, that it must have filled its Roche lobe prior to the second SN. Based on the small proper motion measurement and close proximity to the galactic disk, Piran & Shaviv (2005) conclude that the He star progenitor of J0737−3039B was probably <2.1 M_☉, and the system received a SN kick ≲ 80 km s⁻¹ at birth. This possibility, although more complex in its origin than claimed by Piran & Shaviv (2005), was confirmed by other studies (Willems et al. 2006; Stairs et al. 2006; van den Heuvel 2007; Wong et al. 2010; Ferdman et al. 2013; Dall'Osso et al. 2014). Limits on the He star progenitor mass for this system preclude the second born NS from having ever formed an Fe core (Nomoto 1987); prior to core-collapse, it was only massive enough to form an ONe core. This constraint combined with its very low natal kick requires that this system formed in an ECS. Recently, a population synthesis study by Dewi (2010) showed that the He star progenitor of J0737−3039B could not have overflowed its Roche lobe in a second CE event; instead the mass transfer phase by the He progenitor must have been stable. Thus, we can conclude that J0737−3039 was formed through Channel III.

B1534+12. This system, too, has been studied quite extensively (Willems et al. 2004; Thorsett et al. 2005; Stairs et al. 2006; Wong et al. 2010). All groups agree that the only possible way to explain all the system parameters was if the secondary had gone through a phase of stable mass transfer as a He star. However, this system differs from J0737+3039 in two ways. First, the second NS has a mass of 1.35 M_☉, consistent with the expected remnant of an Fe core-collapse SN. Second, the natal kick magnitude for this second NS is constrained to be high, on order of several hundred km s⁻¹. Based on our assumptions about ECS, we place B1534+12 within Channel II.

B1913+16. Although many authors have commented on the restrictions of the kick velocity and direction required to form B1913+16, the evolutionary constraints on the system's formation channel are less strong (Fryer & Kalogera 1997; Willems et al. 2004; Ihm et al. 2006; Wong et al. 2010). According to these studies, the kick velocity forming the second NS must have been large (≳ 300 km s⁻¹), and the mass of the secondary NS is the expected result of an Fe core-collapse SN, so an ECS is ruled out for the formation of the secondary NS. However, since the mass transfer history of this system has not yet been constrained, we can only restrict B1913+16 to have evolved through either Channel I or Channel II.

To quantitatively analyze the comparative goodness of fit for each model, we adapt the Bayesian analysis used in Ihm et al. (2006) to fit for two independent parameters, orbital period and eccentricity. We begin by applying Bayes's Theorem to our results:

$$\begin{equation} {P}(M|D)={{P}(D|M){P}(M) \over {P}(D)}, \end{equation} \tag{ 1 }$$

where P(M|D) is the probability of the model being correct given the data, P(D|M) is the probability of the data given the model, P(M) is the prior probability of the model, and P(D) is a normalizing constant that is independent of the model, M. Our analysis, however, includes prior evolutionary constraints, so Bayes's Theorem becomes

$$\begin{equation} {P}(M|D,E) = {{P}(D|M,E){P}(E|M){P}(M)\over {P}(D,E)}, \end{equation} \tag{ 2 }$$

where P(M|D, E) is the value which we are calculating: the probability that our model is correct, given the data. P(D|M, E) is the probability of the observed values given our model. P(E|M) is the probability that, given our model, the evolutionary constraints for the data are satisfied. P(M) is our prior, the a priori probability that each particular model is correct. We give each model the same prior probability. Let

$$\begin{equation} C={{P}(M) \over {P}(D,E)} \end{equation} \tag{ 3 }$$

yielding:

$$\begin{equation} {P}(M|D,E)=C{P}(D|M,E){P}(E|M). \end{equation} \tag{ 4 }$$

Because we have given each model the same prior probability, C is independent of model. Let Λ(D)_evol be

$$\begin{equation} \Lambda (D)_{\rm evol} \equiv {P}(D|M,E){P}(E|M). \end{equation} \tag{ 5 }$$

Since C is independent of model, the values of Λ(D)_evol give the relative probabilities of each model. Now, we substitute orbital period, P, and eccentricity, e, for our data in our equation for Λ(D)_evol. Due to the independence of each observed DNS, the probability of the data set given the model is equal to the product of the probabilities of each independent system:

$$\begin{equation} \Lambda (D)_{\rm evol}=\prod _i {P}(e_i,P_i|M,E_i){P}(E_i|M). \end{equation} \tag{ 6 }$$

It is necessary to include a subscript with E because there are different evolutionary constraints on each observed system. P(E_i|M) can be easily determined as the fraction of systems that go into each evolutionary channel.

The errors for the values of eccentricity and orbital period are extremely small and can be ignored. Therefore, calculating P(e, P|M, E_i) reduces to finding the probability density of our model distributions at each point in the 2D space corresponding to the observed DNS. Because our data space is 2D, and each model carries no additional parameters, we do not need to use MCMC-type methods in the computation of the model posterior; instead we calculate P(e, P|M, E_i) by binning the simulated systems into two dimensions for each evolutionary channel. We use a bin size of Δe = 0.033 in eccentricity and Δlog(P) = 0.167 in orbital period. Varying the bin size gives a good handle on the inherent error of our statistical method. Using Equation (5), the value we are calculating, Λ(D)_evol, is equal to the product of the probability density of each bin that corresponds to the observed systems. We find that, although our models have values of Λ(D)_evol that vary by several orders of magnitude, Λ(D)_evol changes within a single model by at most a factor of a few as the bin size varies over reasonable ranges, so our model rankings are insensitive to our choice of binning.

It is useful to determine the dependencies of our results on our use of evolutionary constraints. To do this, we marginalize over the nuisance parameter of the evolutionary channels. The analysis is performed similarly, the only exception being our entire model distribution is binned in two dimensions without evolutionary constraints. In this case, Λ(D) can be defined as

$$\begin{equation} \Lambda (D) = \prod _i {P}(e_i,P_i|M). \end{equation} \tag{ 7 }$$

The values for Λ(D) (hereafter Λ) will, in general, differ from the values we find for Λ(D)_evol (hereafter Λ_evol).

The models in Figures 2 and 3 demonstrate the influence of the energy efficiency in CE evolution, SN kicks, as well as hypercritical NS accretion and survival through the CE phase when the donor star is in the HG.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As expected, tighter orbits are produced by models with less efficient CE phases, those with smaller α_CEλ values. This is reflected in Models 2, 3, and 4 with α_CEλ = 0.25, 0.3, and 1.0, respectively in Figure 2. This figure further demonstrates the effect on branching ratios between the models; smaller CE efficiencies lead to proportionately fewer systems becoming DNSs through Channel III. The interplay of these two effects demonstrate the difficulty in finding models that reproduce J0737−3039. Models 4 and 5, which also have difficulty of producing systems similar to J0737−3039, nonetheless remain viable models due to their increased ability to reproduce other DNSs, particularly J1518+4904 and J1829+2456.

Smaller NS natal kick velocities result preferably in less eccentric DNS orbits. Model 8 with a kick velocity dispersion of 50 km s⁻¹ in Figure 3 shows how the highest concentrations of DNSs tend to have eccentricities smaller than 0.5. Although these low kick models better reproduce the small eccentricity DNSs, they have difficulty reproducing B1913+16. Low kick models alter the branching ratios, creating more DNS through Channel I. Model 9 in Figure 3, which applies kicks with a 150 km s⁻¹ dispersion, shows an intermediate between Model 8 and the standard model.

Model 5 in Figure 2 demonstrates the shift toward slightly longer orbital periods when NSs are allowed to accrete hypercritically in a CE. Such accretion increases the NS mass, reduces the envelope mass to be ejected and therefore requires a smaller degree of orbital contraction for envelope ejection. On the other hand, allowing HG stars to survive the CE phase, if there is enough orbital energy initially, has an almost negligible effect as seen by Model 6 in Figure 3. There is a slightly higher density of systems at the shortest orbital periods in Channel III. Model 7, which both allows NSs to accrete hypercritically as well as allows HG stars survive a CE, shows the combination of these two effects, as expected.

In their study of double compact object merger rates, Dominik et al. (2012) find that only relatively close binaries will overfill their Roche lobe while the donor is on the HG. If these donors are allowed to survive a CE, they will produce tight binaries that merge relatively quickly due to gravitational wave radiation (as evidenced by their short calculated delay time distribution), resulting in an increased DNS merger rate. Since these systems merge soon after their formation, they are unlikely to be observed as a DNS. Therefore, while this evolutionary channel may significantly contribute to the DNS merger rate, its effect on the distribution of observed DNSs in P − e space is negligible.

Using a range for M_{He, deg} of 1.83–2.5 M_☉ populates Channel III, while the number of systems formed through Channel III drops precipitously for a slightly lower range of 1.83–2.25 M_☉. We further find that while we test several different ranges, we find that only those models with a range for M_{He, deg} of 1.83–2.5 M_☉ and 2.0–2.5 M_☉ produce DNSs through Channel III in reasonable numbers. If the He star mass range is not high enough, either the partially degenerate CO core will not become massive enough (>1.08 M_☉) to begin stable C burning, or the resulting ONe core will not reach the Chandrasekhar limit (1.38 M_☉) and leave an ONe WD remnant.

Constraints on M_{He, deg} were independently proposed by Linden et al. (2009) who analyzed the population of high mass X-ray binaries in the Small Magellanic Cloud. They found that if X-ray binaries were not formed through an ECS, the numbers of observed systems with ages 20–60 Myr could not be explained. However, when using a range for M_{He, deg} that allowed X-ray binaries to be formed through ECS, X-ray binaries were formed with the correct age distribution. Although they did not determine the best range for M_{He, deg}, their fiducial range of 1.83–2.25 M_☉ compares to our preferred ranges for M_{He, deg} of 1.83–2.5 M_☉ and 2.0–2.5 M_☉.

Following the analysis of Section 3.3, we evaluate each of our 155 models by calculating the two likelihoods, Λ and Λ_evol, without and with the evolutionary constraints, respectively; the values are shown in Figure 4 in decreasing order across model number. The ordering is performed separately in each panel. When evolutionary constraints are included in the analysis, the number of viable models drops from >100 to 70. We select these 70 models as their Λ_evol value is within three orders of magnitude of the highest likelihood model. Our remaining discussion will be confined to these "best-fit" models indicated by the dotted line in the bottom panel of Figure 4. Table 2 shows the input parameters of the top 70 models, including the reference model.

Zoom In Zoom Out Reset image size

Within the top 70 models, there is little preference either for or against models allowing hypercritical accretion in the CE or allowing HG stars to survive a CE. Furthermore, the secondary parameters have little effect on the resulting distributions. These include changing the mass beyond which a helium core is fully convective, using an angular momentum prescription to calculate the CE evolution, altering the efficiency of mass accretion during stable mass transfer, altering the wind mass loss prescription and decreasing the maximum allowed NS mass.

Two different ranges for M_{He, deg} dominate the top models in Table 2: 1.83–2.5 M_☉ and 2.0–2.5 M_☉. This is due to the difficulty for other ranges to create DNSs through Channel III, a necessity for forming J0737−3039.

There is an interplay between α_CEλ and the natal kick velocities applied in our models: low kick models push the distributions toward smaller eccentricities, while high CE efficiency models push the distributions toward longer orbital periods. Combined, these models make it difficult to form J0737−3039 through Channel III. Therefore, when evolutionary constraints are applied to the models, the models most often eliminated are those with α_CEλ = 1 and a low kick velocity. This is consistent with independent but similar constraints derived by Linden et al. (2009) who argue that large Fe core-collapse kicks and small ECS kicks are required to explain the population of high mass X-ray binaries within the Small Magellanic Cloud. Nevertheless, there are viable models with all three kick velocities that we test. We note that in models with a kick dispersion velocity of 300 km s⁻¹, the bulk of DNSs are formed with e > 0.4, while only two of the eight known DNSs have such high eccentricities. We show in Section 5.2 that this is unlikely to be due to a bias against detecting high eccentricity DNSs.

Figure 5 shows the branching ratios of the top models, ordered by decreasing Λ_evol. The majority of DNSs in these top models are formed through channel II shown as a green line. These systems are characterized by partially recycled primary NSs, and secondary NSs with masses of ∼1.35 M_☉, consistent with many of the DNSs in Table 1. The red line in Figure 5 shows that in these models, DNSs are formed in appreciable numbers through Channel III.

Zoom In Zoom Out Reset image size

One motivation for the present study is to examine whether we can identify the statistically favored evolutionary models for the DNS systems, for which, at present, it is impossible to identify their evolutionary history based on their individual studies. With this in mind, we re-assess the branching ratios for the best-fit Λ_evol models (which have been "calibrated" against the systems with known evolutionary channels) in direct connection with the observed systems. For each observed system and for each best-fit model, we select the sub-population of model systems with eccentricities and orbital periods in the same 2D bin as is used to calculate our Λ_evol values in Section 3.3. The evolutionary channel branching ratios are calculated as the proportion of systems falling within this bin that go through each evolutionary channel. Varying the bin size affects the evolutionary channel branching ratios at the 10%–15% level.

The results for all 70 viable models are shown in Figure 6. For systems such as B1913+16 and J1518+4904, one channel is clearly favored statistically, although this is not necessary the channel that actually formed each system. A few of these models are peculiar. Models at 31 and 48 are models that calculate CE evolution based on the conservation of angular momentum not energy. Although our statistical analysis gives each model an equal prior, modern discussions generally disfavor this prescription (Webbink 2008; Woods et al. 2011; Ivanova et al. 2013).

Zoom In Zoom Out Reset image size

J0737−3039. This is the prime example illustrating that the statistically favored channel is not necessarily the true channel. In this case, the former is Channel II, while the latter is Channel III, based on the analysis of all its observed characteristics (as discussed in Section 3.2). In Figure 6, the red line representing Channel III shows that in all models, no more than half of all DNSs are formed through this pathway.

B1534+12. The requirement of B1534+12 to form through Channel II is in agreement with the branching ratio results shown in Figure 6; in this case, the statistically favored channel is also the true channel. As discussed in Section 3.2 the mass measurements and NS kick constraints indicate that B1534+12 could not have been formed in an ECS, precluding formation through Channel III, the other likely evolutionary channel in Figure 6.

B1913+16. Prior evolutionary constraints restrict B1913+16 to be created through either Channel I or II. In Figure 6 we see that, of these two, Channel II is greatly favored from a statistical point of view. The models at 31 and 48 correspond to models that calculate CE evolution based on an angular momentum conservation description. Hence we can say that B1913+16 most likely went through a phase of stable mass transfer prior to the second SN event. This can be attributed to both Channel II having a much higher branching ratio as well as the systems going through Channel I being formed at larger orbital periods.

J1756−2251. Recent observations by Ferdman et al. (2014) indicate the pulsar companion in J1756−2251 has a mass of 1.23 M_☉. The low mass combined with their newly measured tangential velocity of ∼20 km s⁻¹ suggests that the pulsar companion may have been born in an ECS. If that is the case, J1756−2251 must have been formed through Channel III. The branching ratio for Channel III (red in Figure 6) varies between 10% and 60% for J1756−2251, indicating this possibility. If we add the constraint requiring J1756−2251 to be formed through Channel III, the resulting Λ_evol values do not change significantly.

J1906+0746. J1906+0746 is unique in that the observed pulsar is the second-born NS in the system. A mass measurement for the pulsar indicates it may have been born in an ECS, constraining the system to Channel III. Figure 6 indicates that this system almost certainly went through either Channel II or III. Through future analysis of the galactic position and velocity of J1906+0746, it may be possible to determine if it went through a phase of stable mass transfer, and therefore Channels II or III. However, with the available information we cannot differentiate between the two channels, as the branching ratios in Figure 6 indicate both are possible.

J1811−1736. The branching ratios in Figure 6 indicate that J1811−1736 likely evolved through Channel I or II. Its relatively large orbital period means formation through Channel III is rare for most models. Models 7, 8, and 9 in Figure 3 show this is because of the difficulty of forming a DNS with large orbital periods and eccentricities through the low velocity kick applied in an ECS. The models that allow for the formation of J1811−1736 through Channel III tend to be those with higher kick velocities, and hence larger ECS kicks. Based on their measurements of the relativistic periastron advance and the low derived system mass, Corongiu et al. (2007) suggest that J1811−1736 was born with a low velocity kick. Although uncertainties in the current mass measurements for the system are too large to determine evolutionary constraints, Corongiu et al. (2007) calculate that a second post-Keplerian parameter may be measureable in the near future. If future observations indicate this system was formed through Channel III, it could provide an important constraint on ECS kicks.

J1518+4904. The ∼8.6 day orbital period and eccentricity of 0.249 for J1518+4904 make its formation difficult. The distributions in Figures 2 and 3 show that most DNSs are formed at either smaller orbital periods or larger eccentricities. Although a few models form J1518+4904 through Channel II or III, Figure 6 indicates it was most likely formed through Channel I. Observational errors on the mass estimates are too large to indicate whether the companion was born in an ECS event or not.

J1829+2456. Depending on the individual model, Figure 6 shows that J1829+2456 formed through any of the three evolutionary channels. Poor mass constraints do not allow any further indication. Future observations may provide further constraints on its evolution.

While there is little observational evidence to distinguish between different evolutionary channels, one such piece of evidence could be the degree to which the first-born NS has been recycled. Ivanova (2011) has argued that the mass losing star in a CE will go through a phase of thermal readjustment upon envelope expulsion. During this phase an NS companion can be mildly recycled. Investigation of Table 1 shows that there may be in fact a bifurcation spin periods in the primary NS: those with P_s ≲ 50 ms and those with P_s ≳ 100 ms. If this argument is accurate, then J1811−1736 went through Channel I, avoiding any stable mass transfer after the CE phase, while the other systems with shorter spin periods went through Channel II or III. J1906+0746 is an exception because we do not know the spin period of the primary NS.

Comparison with Figure 6 shows that this idea is consistent with most of the results of our analysis here. In roughly half of our viable models, Channel I is preferred for J1811−1736. With a few exceptions, other systems prefer Channel II (in agreement with their shorter spin periods) except for J1518+4904 which our results here indicate was likely formed through Channel I. Either this system, too, formed through Channel II, or the primary NS in J1518+4904 accreted enough material in a CE to recycle the primary to its current spin period of 40.9 ms.

Our statistical analysis discussed in Section 3.3 relies on the assumption that all DNSs have an equal detection likelihood. For example, a strong bias against detection of DNSs with large eccentricities means that high kick velocity models, which produce on average more eccentric DNSs, will not be as disfavored as the Λ_evol values indicate. In principle our DNS distributions could be convolved with any observational bias for a more accurate comparison to the observed DNSs.

One such bias deals with the observability of pulsars in DNSs which is a function of a number of intrinsic parameters such as the beaming fraction, magnetic field strength, age, and spin period as well as particular parameters such as the distance and direction from the Sun. However, since we ignore non-interacting DNSs, our model assumes that every DNS contains a partially recycled pulsar, each with an equal probability of being observed. Absent any deeper understanding of the connection between a pulsar's characteristics and its prior evolution, we take this to be a reasonable approximation.

An additional bias could be caused by Doppler smearing of pulsars in short orbital periods. Due to the high accelerations in such systems over the course of a single observation, close binaries suffer a decreased detection efficiency. Bagchi et al. (2013) analyzed the detectability of binary pulsars in surveys, finding that modern acceleration searches (e.g., PRESTO) were greater than 80% efficient at detecting pulsars in DNSs for typical parameters. Interestingly, they found systems with the shortest orbital periods and lowest eccentricities suffered the strongest biases, which still have a detection efficiency greater than 50%.

A potentially more pernicious bias could be introduced by the limited sky coverage of pulsar surveys, typically close to the Galactic plane. Binaries that gain large systemic velocities due to the NS natal kicks may therefore be missed. Such DNSs would likely have significant eccentricity and a large orbital period, possibly helping to explain the relative dearth of such systems despite the predictions by our simulations. We first investigate the systemic velocities of our resulting DNS populations.

The distribution of systemic velocities, normalized to the highest bin, for Models 8 (σ_kick = 50 km s⁻¹), 9 (σ_kick = 150 km s⁻¹), and 1 (σ_kick = 300 km s⁻¹) are shown in Figure 7. We expect these three models to be representative of the different kick velocity models of our entire set. Calculated from the equations in Kalogera (1996), this systemic velocity arises from both the mass lost from the collapsing object as well as the kick applied to the newly born NS. Model 8 produces DNSs with systemic velocities peaking at 35 km s⁻¹. Model 9 and 1 both produce DNSs with systemic velocities peaking at 75 km s⁻¹, although the higher kick velocity Model 1 has a high velocity tail. Can these systemic velocities create a bias against detection strong enough to significantly affect our quantitative results?

Zoom In Zoom Out Reset image size

To determine this magnitude of any effect, we create a model Milky Way-like galaxy composed of a three-component potential defined in Appendix A of Wex et al. (2000). DNSs are born in a double exponential, axisymmetric disk with a scale length of 2.80 kpc and scale height of 70 pc, and are given an initial velocity corresponding to circular rotation. The motion of the DNS through the Galaxy is calculated using a Runge–Kutta fourth-order scheme. The binary is evolved for 20 Myr (the approximate evolution time of the binary prior to the birth of the second NS) before a systemic velocity is applied. The systemic velocity is applied to the system in a randomly chosen direction. The system is then evolved for a Gyr and the location of the binary recorded throughout this evolution at intervals of a Myr. This location is translated into a galactic longitude and latitude as referenced from the Sun, placed at 8.5 kpc from the center of the Galaxy. We repeat this process, simulating the motion of 1000 binaries throughout the Galaxy, each with a different randomly chosen initial position and systemic velocity direction. We then determine the percentage of the resultant Galactic longitudes and latitudes that fall within a region near the Galactic Plane, using as a fiducial region the Parkes multibeam pulsar survey area (Manchester et al. 2001): |b| < 5° and 260° < l < 50°. Figure 8 shows this percentage over the relevant range of systemic velocities. We find that when no velocity is applied, a DNS has a 93% chance of falling within the Parkes multibeam pulsar survey region. That number falls to 85% for systemic velocities of 100 km s⁻¹ and 69% for systemic velocities of 200 km s⁻¹.

Zoom In Zoom Out Reset image size

We test the strength of this bias on Models 8, 9, and 1 with Maxwellian kick dispersion velocities of 50, 150, and 300 km s⁻¹, respectively. For each simulated system in the three models, we find the detection probability from Figure 8 corresponding to its individual systemic velocity. We then convolve observational biases into our distributions in P − e space by selecting a subset of our original distribution. Each simulated DNS has a chance equal to its calculated detection probability of being included into the subset. In Model 1 (with the highest kick velocities and hence the strongest expected effects of the three models) roughly 80% of all systems are unaffected. The effects of this observational bias are seen in Figure 9; the top panels provide the full distributions of DNSs for our three test models while the bottom panels show the distributions convoluted with the observational bias. The resulting distributions are nearly qualitatively identical. Since our quantitative results are robust to a factor of several orders of magnitude, we can safely ignore this observational bias as potentially altering our conclusions in this study.

Zoom In Zoom Out Reset image size