Formation of Double Neutron Star Systems

In our galaxy, the Milky Way, it is estimated that some ${10}^{8}\mbox{--}{10}^{9}$ NSs are present. As a rough estimate, one can simply assume a constant NS formation rate of 2 per century for 10 Gyr. However, we can only detect a tiny fraction of these—namely, those that are rapidly spinning and possess strong magnetic fields (thereby producing radio pulsars that may beam in our direction), those that are accreting material from a companion star in a close binary (giving rise to the emission of X-rays), or those that are nearby, young, and hot (giving rise to thermal emission). According to the ATNF Pulsar Catalogue¹¹ (Manchester et al. 2005), there are $\sim 2600$ Galactic radio pulsars detected at present (2017 June), out of a total known population of roughly 3000 NSs in the Milky Way. The two Magellanic Clouds host some $\sim 30$ known radio pulsars (Crawford et al. 2001; Ridley et al. 2013), besides an impressive number of $\sim 140$ HMXBs (Coe & Kirk 2015; Antoniou & Zezas 2016), some of which are potential progenitors of DNS systems.

Among the new DNS systems recently discovered, the most notable are PSRs J1930−1852 (Swiggum et al. 2015), J0453+1559 (Martinez et al. 2015), and J1913+1102 (Lazarus et al. 2016). PSR J1930−1852 has an orbital period of 45 days and is by far the widest DNS system known. At the same time, this radio pulsar has the slowest spin period (185 ms) of any recycled pulsar in a DNS system. This system thus provides important data for understanding the formation of wide-orbit DNS systems. PSR J0453+1559 is a DNS system with an asymmetric pair of NS masses (1.56 and $1.17\,{M}_{\odot }$), deviating significantly from the "canonical" value of $\sim 1.35\,{M}_{\odot }$ (Thorsett & Chakrabarty 1999). This system is the first to finally break out of the narrowly limited observed NS mass interval, which was previously found to be $1.23\mbox{--}1.44\,{M}_{\odot }$ for all NSs in DNS systems, thereby alleviating a theoretical conundrum on NS birth masses (e.g., Schwab et al. 2010; Tauris et al. 2011, and references therein). A third new DNS binary is PSR J1913+1102 (Lazarus et al. 2016), which has the largest total mass of any known DNS system, thus potentially yielding another asymmetric-mass system. Its GW merger timescale ($\sim 480\,\mathrm{Myr}$) is on the short side.

To understand the formation of DNS systems from a binary stellar evolution point of view, it is important to consider all data provided by observations and their derived parameters, e.g., NS masses (${M}_{\mathrm{NS}}$), spin periods (P), their time derivatives ($\dot{P}$) and hence surface magnetic fields (B), orbital periods (${P}_{\mathrm{orb}}$), and eccentricities (e), as well as kinematic properties such as systemic velocities (${v}_{\mathrm{sys}}$) and their Galactic locations.

Geometric constraints (e.g., misalignment angles between the spin axis of the recycled radio pulsar and the orbital angular momentum vector) are also important for probing the SN explosion. Direct relations between DNS systems and potential SN remnants (none detected so far) would be helpful to constrain their ages. In addition, any correlations between, for example, ${P}_{\mathrm{orb}}$, P, and eccentricity are crucial for understanding the recycling process and needs to be investigated throughly—one of the major aims of this paper.

In the vast majority of DNS systems, only the (mildly) recycled radio pulsar is seen; the reason being that the mass accretion process decreases the B-field strength of the NS significantly (Bhattacharya 2002), whereby the loss rate of rotational energy becomes much lower due to weaker braking torques from pulsar winds, magnetic dipole radiation, and plasma currents in the magnetosphere, thus making the NS observable as a radio pulsar on a much longer timescale (Srinivasan & van den Heuvel 1982). This is particularly so because the B-fields of recycled pulsars apparently do not decay further (Kulkarni 1986; van den Heuvel et al. 1986). For these recycled DNS pulsars, the values of P ($23\mbox{--}185\,\mathrm{ms}$) and their derivatives ($\dot{P}=2.2\times {10}^{-20}\mbox{--}1.8\times {10}^{-17}\,{\rm{s}}\,{{\rm{s}}}^{-1}$) are direct measures of the efficiency of the recycling process. The corresponding derived surface dipole B-fields are approximately in the range $0.3\mbox{--}18\times {10}^{9}\,{\rm{G}}$. Table 1 provides the observed ranges of key parameters of DNS systems.

Table 1.  Observed Ranges of Key Properties of DNS Systems

Properties of Recycled (Old) NSs:
Spin period, P	$23\mbox{--}185\,\mathrm{ms}$
Period derivative, $\dot{P}$	$(0.027\mbox{--}18)\times {10}^{-18}\,{\rm{s}}\,{{\rm{s}}}^{-1}$
Surface dipole B-field, B	$(0.29\mbox{--}18)\times {10}^{9}\,{\rm{G}}$
Mass, ${M}_{\mathrm{NS},1}$	1.32–1.56 M⊙a
Properties of Young NSs:
Spin period, P	$144\mbox{--}2773\,\mathrm{ms}$
Period derivative, $\dot{P}$	$(0.89\mbox{--}20)\times {10}^{-15}\,{\rm{s}}\,{{\rm{s}}}^{-1}$
Surface dipole B-field, B	$(2.7\mbox{--}5.3)\times {10}^{11}\,{\rm{G}}$
Mass, ${M}_{\mathrm{NS},2}$	$1.17\mbox{--}1.39\,{M}_{\odot }$
Orbital Properties:
Orbital period, ${P}_{\mathrm{orb}}$	$0.10\mbox{--}45\,\mathrm{days}$
Eccentricity, e	$0.085\mbox{--}0.83$
Merger time, ${\tau }_{\mathrm{gwr}}$	$86\,\mathrm{Myr}\to \infty$
Systemic velocity, ${v}_{\mathrm{sys}}$	$25\mbox{--}240\,\mathrm{km}\,{{\rm{s}}}^{-1}$

Note. Data taken from the ATNF Pulsar Catalogue (Manchester et al. 2005)—see Table 2 for further details. Only DNS systems in the Galactic disk are listed. The systemic recoil velocity, ${v}_{\mathrm{sys}}={v}^{\mathrm{LSR}}$, is quoted with respect to the local standard of rest (Section 2.2).

^a1.32 M_⊙ Mark an upper limit to the lowest mass of the first-born NS.

Download table as:  ASCIITypeset image

The observed orbital periods now span between ${P}_{\mathrm{orb}}=0.10\mbox{--}45\,\mathrm{days}$. As we shall discuss in this paper, selection effects make it somewhat difficult to detect radio pulsars in DNS systems at the lower end of this range, as well as second-born pulsars with slow spin periods of a few seconds. The measured eccentricities (e = 0.08–0.83) and estimated systemic velocities, in the local standard of rest (LSR) reference frame (${v}_{\mathrm{sys}}={v}^{\mathrm{LSR}}\approx 25\mbox{--}240\,\mathrm{km}\,{{\rm{s}}}^{-1}$), are important for constraining the SN physics of the second-born NS and reflect the combined effects of the amount of mass ejected and the kick velocity imparted onto the newborn NS.

Previous theoretical works on the physics of DNS formation includes (here disregarding general population synthesis studies) Bisnovatyi-Kogan & Komberg (1974), Wheeler et al. (1974), Flannery & van den Heuvel (1975), Srinivasan & van den Heuvel (1982), van den Heuvel (1994), Ivanova et al. (2003), Dewi & Pols (2003), Podsiadlowski et al. (2004), van den Heuvel (2004), and Dewi et al. (2005). From these papers, a standard scenario¹² has emerged (e.g., Bhattacharya & van den Heuvel 1991; Tauris & van den Heuvel 2006), which we now summarize in more detail.

In Figure 1, we show an illustration of the formation of a DNS system. The initial system contains a pair of OB-stars that are massive enough¹³ to terminate their lives in a core-collapse SN (CCSN). To enable the formation of a tight DNS system in the end, the two stars must initially be in a binary system close enough to ensure interactions via either stable or unstable mass transfer. If the binary system remains bound after the first SN explosion (which is of Type Ib/c; Yoon et al. 2010), the system eventually becomes observable as a HMXB. Before this stage, the system may also be detectable as a radio pulsar orbiting an OB-star, e.g., as in PSRs B1259−63 (Johnston et al. 1992) and J0045−7319 (Kaspi et al. 1994). When the secondary star expands and initiates full-blown Roche-lobe overflow (RLO) during the HMXB stage, the system eventually becomes dynamically unstable. For wide systems, where the donor star has a deep convective envelope at the onset of mass transfer (i.e., during the so-called Case B RLO, following the termination of core hydrogen burning), the timescale on which the system becomes dynamically unstable might be as short as a few 100 yr (Savonije 1978). This leads to the formation of a CE (Paczyński 1976), where the dynamical friction of the motion of the NS inside the giant star's envelope often causes extreme loss of orbital angular momentum and (in some cases) ejection of the hydrogen-rich envelope. If the binary system survives the CE phase, it consists of a NS orbiting a helium star (the naked core of the former giant star). Depending on the orbital separation and the mass of the helium star, an additional phase of mass transfer (Case BB RLO; Habets 1986; Tauris et al. 2015) may be initiated. This stage of mass transfer is important since it enables a relatively long phase of accretion onto the NS, whereby the NS is recycled, and it allows for extreme stripping of the helium star prior to its explosion (as a so-called ultra-stripped SN; Tauris et al. 2013, 2015; Suwa et al. 2015; Moriya et al. 2017). Whether or not the system survives the second SN depends on the orbital separation and the kick imparted onto the newborn NS (Flannery & van den Heuvel 1975; Hills 1983; Tauris & Takens 1998). As we shall argue in this paper, we expect most systems to survive the second SN explosion. If the post-SN orbital period is short enough (and especially if the eccentricity is large), the DNS system will eventually merge due to GW radiation and produce a strong high-frequency GW signal and possibly a short GRB (e.g., Eichler et al. 1989). The final remnant is most likely a BH, although, depending on the EoS, a NS (or, at least, a metastable NS) may be left behind instead (Vietri & Stella 1998).

Zoom In Zoom Out Reset image size

Aside from the pre-HMXB evolution, which is discussed in Section 3.1, the most important and uncertain aspects of our current understanding of DNS formation are related to

• (i)  
  CE evolution and spiral-in of the NS,
• (ii)  
  momentum kicks imparted onto newborn NSs, and
• (iii)  
  the mass distribution of NSs.

In the following subsections, we provide more details on each of these three aspects.

1.5.1. CE Evolution

1.5.2. NS Kicks

1.5.3. Mass Distribution of NSs

Following this introduction, we present the latest observational data on DNS systems in Section 2, and discuss their observational selection effects and biases, with a special focus on DNS masses and the DNS systemic velocities. In Section 3, we discuss the binary evolution from the ZAMS stage until the HMXB stage. We derive in detail in Section 4 the amount of mass accreted by the first-born NS during its evolution in a binary system, with the aim of inferring the birth masses of these first-born NSs in DNS systems. This allows us to directly compare our results with the masses of second-born NSs and also NS masses in other systems such as HMXBs. In Section 5, we analyze DNS correlations between $({P}_{\mathrm{orb}},P)$ and $({P}_{\mathrm{orb}},e)$, and derive theoretical correlations, based on stellar evolution models and accretion physics, which we compare with observational data. We argue that the suggested $(P,e)$ correlation is a natural consequence of the two aforementioned correlations, and simulate the spread in all three correlations as a result of a kick velocity added to the second-born NS. We analyze and discuss the dynamical consequences of asymmetric SNe (kicks) during SN explosions in Section 6, and also investigate a possible connection between the final stellar structure, kick magnitude, and NS mass. In Section 7, we discuss the mapping of observed DNS systems to simulated DNS systems (for which we compute properties right after the second SN), in terms of parameter evolution over time and observational selection effects. Monte Carlo simulations for all known DNS systems are presented in Section 8 in order to derive their pre-SN properties and place interesting constraints for some of these systems, which are useful to investigate the nature of the second SN explosion—one of the prime aims of this paper. The details of the simulations of each individual DNS system are given in the Appendix. The ramifications of these SN simulations are discussed further in Section 9. Finally, in Section 10, we discuss DNS mergers (with applications to LIGO/VIRGO and observations of short GRBs) in the context of this work. Our conclusions are summarized in Section 11.

Before engaging in an analysis of DNS systems, we must first assess any selection effects associated with discovering DNS systems. Most of the recycled pulsars in DNS systems have spin periods of a few tens of milliseconds. This implies that they are near the optimal sensitivity of current surveys for dispersion measure values (DMs) of the order of $325\,{\mathrm{cm}}^{-3}\,\mathrm{pc}$ and lower (see, e.g., Lazarus et al. 2015). This is not the case for slow, non-recycled radio pulsars in DNS systems like PSR J0737−3039B, for which Lazarus et al. (2015) determined a significant DM-dependent loss of sensitivity caused by gain and system temperature fluctuations in radio receivers, due mostly to radio frequency interference. This is particularly dramatic at low DMs where the loss of sensitivity can be one order of magnitude.

The main selection effect against the detection of a recycled pulsar in a tight DNS system is caused by the orbital acceleration of the pulsar (this does not affect the detection of slow pulsars). The fast-changing Doppler shift then means that the spin period is perceived (on Earth) to change significantly within the time of a single observation. Unless corrected, this effect smears the pulsar signal in the Fourier domain, where the searches are performed; particularly for the higher harmonics, it becomes increasingly severe for cases where the observation length (${T}_{\mathrm{obs}}$) is more than a few percent of the orbital period (${P}_{\mathrm{orb}}$). Several strategies have been adopted to cope with this. These include mostly "acceleration searches" (Anderson et al. 1990; Camilo et al. 2000; Ransom et al. 2002), which assume a constant acceleration during the observation. However, these surveys are computationally expensive. While they increase the number of trials necessary for detecting a binary pulsar, they necessarily cause a reduction in sensitivity even if the acceleration is nearly constant. The reason for this is that the much larger number of trials increases the number of candidates that originate from the Gaussian noise in the data. As a result, the signal-to-noise ratio threshold for a detection in an acceleration search is generally larger than in a survey with no acceleration. Furthermore, these surveys start losing sensitivity when ${T}_{\mathrm{obs}}/{P}_{\mathrm{orb}}\,\geqslant \,0.1$ (this number becomes smaller for larger accelerations), i.e., when the assumption of a constant acceleration no longer holds (Ng et al. 2015). Thus, systems with large eccentricities can also be detrimental. Acceleration searches are now widely used in blind surveys (e.g., Lazarus et al. 2015), because computing improvements make them feasible as demonstrated by their successes in discovering new DNS systems.

Several strategies have been pursued to search for pulsars in even tighter systems; these are too numerous to describe here (see e.g., Eatough et al. 2013, and references therein). They are not yet widely used in large-scale blind surveys because of their computational cost.

To summarize, there is generally a bias against the detection of pulsars in DNS systems, particularly in very tight systems, but this depends crucially on the length of the observations in each survey (the shorter the integration, the smaller is the bias) and the particular search algorithms being used in each survey.

When comparing the outcome of our SN simulations (see Section 8) with observational data, it is necessary to calculate the systemic velocities of the DNS systems relative to their local standard of rest (LSR) in the Galaxy.

From pulsar timing, one often gets an accurate value for the proper motion of the DNS system on the sky. Combined with an estimation of the pulsar distance, a transverse velocity, ${v}_{T}^{\mathrm{SSB}}$, with respect to the solar system barycenter (SSB) can be found. Using a Monte Carlo (MC) simulation, we account for the uncertainties in distance and proper motion. Concerning the pulsar distance, in some cases the only available estimation comes from the DM distance, which is based on a model for the free electron density in our Galaxy (e.g., Cordes & Lazio 2002; Yao et al. 2016). In these cases, we assume in our MC simulation a 20% uncertainty, and use this uncertainty to be 1σ in a Gaussian distribution for the DNS distance.

For DNS systems, one does not have a measurement of the radial velocity, unlike the case for some nearby pulsar+WD systems (see, e.g., Callanan et al. 1998; Antoniadis et al. 2012) for which phase-resolved spectroscopy of the optically bright WD is possible. Therefore, one does not have a 3D velocity vector of the DNS system. At this stage, we make the assumption that the spatial orientation of the velocity vector has no preferred direction with respect to the SSB. Adding a random orientation to our MC simulation, one then obtains a probability distribution for the radial component, ${v}_{R}^{\mathrm{SSB}}$ (based on the calculated ${v}_{T}^{\mathrm{SSB}}$): ${v}_{R}^{\mathrm{SSB}}\,={v}_{T}^{\mathrm{SSB}}\cot \chi$, where χ denotes the angle between the direction to the pulsar and the orientation of the 3D velocity vector.

Finally, to convert the SSB velocity vector to the LSR, we have used the Galactic model of McMillan (2017). This conversion is done for every single MC realization in our simulation. Consequently, we get a probability distribution for the magnitude of the systemic velocity with respect to the LSR, ${v}^{\mathrm{LSR}}$, which we then compare with our SN simulations later in this paper.

The results for the systemic velocities, for those DNS systems near the Galactic plane with well-measured proper motion, are given in Table 3. In the last four columns, we state the median value, the $1\sigma$ and $2\sigma$ ranges, as well as the 90% confidence limit for the upper limit of the 3D systemic velocity, ${v}^{\mathrm{LSR},\max }$, based on our MC simulations (strictly speaking, since the ${v}^{\mathrm{LSR}}$ distributions are non-Gaussian, these $1\sigma$ and $2\sigma$ values refer to the 68% and 95% confidence intervals, respectively). We notice that we make the assumption that the LSR is a reasonable representation of the birth reference frame of each DNS system. That might not necessarily be the case, as for instance discussed for the Hulse–Taylor pulsar in Wex et al. (2000). In cases where no precise pulsar distance is known from parallax measurements or orbital decay, we have applied the DM distances based on the NE2001 model (Cordes & Lazio 2002). We adopt this model as the default¹⁴ in our analysis of DNS kinematics discussed further in this work.

Table 3.  Estimated 3D Systemic LSR Velocities for Galactic Disk DNS Systems with Proper Motion Measurements

	l	b	Distance	μα	μδ	${v}_{T}^{\mathrm{SSB}}$	${v}_{\mathrm{median}}^{\mathrm{LSR}}$	${v}_{1\sigma }^{\mathrm{LSR}}$	${v}_{2\sigma }^{\mathrm{LSR}}$	${v}_{90 \% \,\mathrm{CL}}^{\mathrm{LSR},\ \max }$
Radio Pulsar	(deg)	(deg)	(kpc)	(mas yr−1)	(mas yr−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)
J0453+1559	184.12	−17.14	1.07	−5.5	−6.0	41 ± 8	52	36−85	26−197	102
J0737−3039A	245.24	−4.50	1.15 (a)	−3.8	2.1	25 ± 4	33	18−57	10−124	68
J1518+4904	80.81	54.28	0.63	−0.7	−8.5	26 ± 5	28	20−49	15−123	61
B1534+12	19.85	48.34	1.05 (b)	1.5	−25.3	127 ± 1	136	118−224	115−512	276
J1756−2251	6.50	0.95	0.73 (c)	−2.4	5.5	21 ± 7	41	21−71	11−123	82
B1913+16	49.97	2.12	9.80 (d)	−1.2	−0.8	69 ± 22	240	174−273	83−401	290
J0453+1559	184.12	−17.14	0.52	−5.5	−6.0	20 ± 4	27	20−43	15−99	52
J1518+4904	80.81	54.28	0.96	−0.7	−8.5	39 ± 8	39	28−70	21−184	89
B1913+16	49.97	2.12	5.25	−1.2	−0.8	37 ± 7	132	98−168	67−229	179

Note. The upper part of the table uses NE2001 for the DM distances to PSRs J0453+1559 and J1518+4904; otherwise, distances are taken from (a) Deller et al. (2009), (b) Fonseca et al. (2014), (c) Ferdman et al. (2014), and (d) Weisberg et al. (2010). The lower part of the table uses Yao et al. (2016) for DM distance determinations. See Section 2.2 for further details.

Download table as:  ASCIITypeset image

The lower part of Table 3 shows our results using the DM distances based on Yao et al. (2016), who suggest that their model is better than NE2001. However, more thorough studies of this new model are needed before any firm conclusions can be drawn. For some DNS systems, like PSR J0737−3039, the results obtained from the different DM distance models are quite similar, whereas in other cases (PSRs J0453+1559 and B1913+16) the value of ${v}^{\mathrm{LSR}}$ can differ by a factor of two. One particular example of a DNS pulsar with a discrepancy is PSR J1756−2251 where Yao et al. (2016) find a DM distance of 2.81 kpc, although an upper value of 1.2 kpc for its distance can be derived (Ferdman et al. 2014) using orbital period decay (${\dot{P}}_{\mathrm{orb}}$).

To further check our results against the sensitivity of the DM distance, we recalculated the values of ${v}^{\mathrm{LSR}}$ for the Yao et al. (2016) DM distances, now assuming in our MC simulation a 50% uncertainty in distance using a flat probability range centered on the DM distance. The resulting ranges of values for ${v}_{1\sigma }^{\mathrm{LSR}}$ and ${v}_{2\sigma }^{\mathrm{LSR}}$ only changed by a few $\mathrm{km}\,{{\rm{s}}}^{-1}$ for PSRs J0453+1559 and J1518+4904, whereas for PSR B1913+16 these interval ranges widened on each end by $\sim 18\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and $\sim 12\,\mathrm{km}\,{{\rm{s}}}^{-1}$, respectively. Thus, these are relatively small changes compared to the values stated in Table 3 and therefore not significant for our NS kick discussions in Section 6.

To estimate the (equatorial) surface B-fields of pulsars, the standard method is to use the classic vacuum magnetic dipole expression (e.g., Lorimer & Kramer 2004): $B\simeq 3.2\,\times {10}^{19}\,{\rm{G}}\,\sqrt{P\dot{P}}$. The incompleteness of the vacuum magnetic dipole model, however, is particularly evident after the discovery of intermittent pulsars by Kramer et al. (2006a) and demonstrates the need for including the plasma term in the spin-down torque. The classical expression does not include the rotational energy loss contribution from the pulsar wind nor the spin-down torque caused by the ${\boldsymbol{j}}\times {\boldsymbol{B}}$ force exerted by the plasma current in the magnetosphere (e.g., Goldreich & Julian 1969; Kramer et al. 2006a; Spitkovsky 2008). Hence, the classic expression does not predict any spin-down torque for an aligned rotator ($\alpha =0^\circ$), which is not observed. A combined model was derived by Spitkovsky (2006; see also Li et al. 2012), and applying his relation between B and α, Tauris et al. (2012) obtained¹⁵

$$\begin{eqnarray}B & = & \sqrt{\displaystyle \frac{{c}^{3}{IP}\dot{P}}{4{\pi }^{2}{R}^{6}}\,\displaystyle \frac{1}{1+{\sin }^{2}\alpha }}\\ & \simeq & 1.3\times {10}^{19}\,{\rm{G}}\quad \sqrt{P\dot{P}}\,{\left(\displaystyle \frac{M}{1.4{M}_{\odot }}\right)}^{3/2}\sqrt{\displaystyle \frac{1}{1+{\sin }^{2}\alpha }}.\end{eqnarray} \tag{ 1 }$$

To estimate the surface B-fields of our DNS pulsars, we use the above expression for an assumed magnetic inclination angle $\alpha =60^\circ$. This yields roughly the following relation: $B\simeq 1.0\,\times {10}^{19}\,{\rm{G}}\,\sqrt{P\dot{P}}$. Figure 2 shows the distribution of DNS radio pulsars in the ($P,\dot{P}$) diagram and their estimated B-fields.

Zoom In Zoom Out Reset image size

Precise measurements of NS masses serve as an important diagnostic tool for understanding NS formation. In the following subsection, we therefore outline in detail how to determine such NS masses from observations.

There are a wide variety of methods available for measuring NS masses (Özel & Freire 2016). For DNS systems specifically, all mass measurements to date result from the measurement of relativistic post-Keplerian (PK) parameters via high-precision radio pulsar timing. In DNS binaries, where classical tidal or magnetic effects are absent, PK parameters describe relativistic corrections to the timing and pulse-structure data. Most notable are effects related to the orbital motion and time dilation of the pulsar and the propagation of the pulsar signals in the curved spacetime of the binary system (Damour & Taylor 1992). Generally, for a given theory of gravity, the PK parameters are (to leading order) functions only of the (well-known) Keplerian orbital parameters and the masses of the two components of the system. Hence, a measurement of two PK parameters determines the masses of the pulsar and companion uniquely. In several DNS systems, more than two PK parameters have been measured, which then also allows for a consistency test of the applied gravity theory. So far, general relativity (GR) has passed all these tests with flying colors (Wex 2014), and has therefore been confirmed to be the correct theory to determine the masses from PK parameters in DNS systems.

The validity of this approach has been particularly well demonstrated for the double pulsar system PSR J0737−3039A/B, where six PK parameters have been measured so far, plus (uniquely in a DNS) the component mass ratio provided by the detection of the second NS as a radio pulsar. Crucially, all 21 pairwise combinations of mass ratio and PK parameters yield consistent (and in most cases very precise) mass values for the two pulsars in this system. These agreements can be reframed as a total of five high-precision tests of GR, with one PK parameter and the mass ratio providing the pulsar masses, and the remaining PK parameters representing one test of GR each (Kramer et al. 2006b; Breton et al. 2008). Aside from confirming the validity of GR, these measurements also confirm that there are no unknown additional classical effects perturbing the timing of the pulsar in a measurable way.

Although the vast majority of MSP+WD systems generally have higher timing precision than DNS systems, they often have poorly measured PK parameters. One of the reasons for this is their orbital eccentricity, which is many orders of magnitude smaller than for most DNS systems. A large eccentricity implies that the system has a well-defined periastron, allowing for the measurement of its rate of advance, $\dot{\omega }$, which immediately yields the sum of the masses of the two components of the system in GR. The orbital eccentricity also allows the measurement of the "Einstein delay," γ (which is a combination of variation in the special relativistic time dilation, due to the changing orbital velocity of the pulsar, and a variation in the gravitational redshift for the pulsar due to a changing distance to the companion). However, the latter effect can be re-absorbed in the definition of the size of the Keplerian orbit. For it to be separately measurable, the orbit has to precess by at least a few degrees—something that would take many decades to millennia for the wider DNS systems to be measured with present radio telescopes.

The above two PK parameters ($\dot{\omega },\gamma$) provide most of the precise NS mass measurements in tight DNS systems. The remaining PK parameters can be observed even for circular orbits. The emission of GWs gives rise to an orbital decay, ${\dot{P}}_{\mathrm{orb}}$. Like γ, this can only be measured for the more compact systems. A large orbital eccentricity magnifies the effect significantly (Peters 1964).

Finally, the Shapiro delay is a relativistic effect that yields two PK parameters: r and s in the Damour & Deruelle (1986) parameterization, or h₃ and ς in the Freire & Wex (2010) parameterization. These allow for a direct determination of the mass of the pulsar's companion ($r={{GM}}_{\mathrm{comp}}/{c}^{2}$) and the orbital inclination i, more precisely $\sin i$. The precision of this method is independent of the orbital period, but it requires a fortunate combination of high orbital inclination and good timing precision.

For the wide-orbit DNS systems, γ and ${\dot{P}}_{\mathrm{orb}}$ are generally not measurable and a non-detection of the Shapiro delay implies that we cannot separate the two NS masses. In such cases, only the total mass of the system is known from $\dot{\omega }$.

In Figure 3, we have plotted the distribution of NS masses in both DNS and NS+WD systems. The data are taken from Özel & Freire (2016) and Table 1 in Antoniadis et al. (2017). Whereas DNS systems descend from HMXBs, NS+WD systems are usually produced in LMXBs. The error bars on NS masses from the precisely measured DNS systems are $\lesssim 0.01\,{M}_{\odot }$, whereas the error bars for the pulsar masses in the NS+WD systems plotted here range between $0.1 \% \mbox{--}10 \%$. Two NS+WD systems (PSRs B2303+46 and J1141−6545) are unique in the sense that here the NS is believed to have formed after the WD (Tauris & Sennels 2000). These NSs are observed with the typical characteristics of young and non-recycled radio pulsars (i.e., large values of P and $\dot{P}$), and further, their orbits are eccentric and not circular as in the vast majority of MSP binaries. For these two systems, we take the NS mass estimates of van Kerkwijk & Kulkarni (1999) and Bhat et al. (2008). It is expected that recycled pulsars will have larger masses than non-recycled pulsars since a certain amount of accreted mass is needed to recycle a pulsar. Indeed, for the DNS systems (Figure 3, top and central panels) it seems that the recycled NSs might have larger masses by an amount of the order $\sim 0.1\,{M}_{\odot }$. However, a crucial question is how much mass these NSs actually accreted, and whether the mass difference between second-born (non-recycled) NSs and recycled NSs is instead mainly caused by differences in the birth masses of NSs resulting from the different conditions in the first and the second SN explosion, respectively. We investigate this question in more detail later in this paper.

Zoom In Zoom Out Reset image size

It has been shown recently that interacting with a binary companion during its lifetime is the rule rather than the exception for a massive star. The majority of them are found in close binaries that will start mass transfer before the first SN occurs (Kiminki & Kobulnicky 2012; Sana et al. 2012; Sota et al. 2014).

If the mass transfer (i.e., RLO) remains stable, it results in the almost complete loss of the hydrogen-rich envelope of the initially more massive star, a significant fraction of which will be accreted by its companion (Podsiadlowski et al. 1992); see Figure 1. The accretion process leads to a spin-up of the mass gainer, which is expected to end up spinning close to critical rotation (Packet 1981; Habets 1986; de Mink et al. 2013) and to a net widening of the orbit. Most models predict that the mass donor will reach the SN stage first. If the system does not break up due to the SN (Section 6), its observational counterparts include the so-called Be-star X-ray binaries (BeHMXBs; Section 3.2), where the Be nature of the core-hydrogen burning component signals its state of near-critical rotation (Townsend et al. 2004). The major uncertainty in our understanding of this evolutionary path concerns the mass and angular momentum loss from the system during the mass-transfer phases, i.e., the amounts, but also the temporal behavior, of mass and angular momentum loss (Langer 2012). Some BeHMXBs may later evolve into supergiant HMXBs (sgHMXBs), similar to those observed in wide orbits (Section 3.2).

In case the binary mass transfer becomes unstable, the system will undergo a CE phase where both components share the hydrogen-rich envelope of the primary. There are two possible outcomes from this (Podsiadlowski et al. 1992). One is that the two stars merge into a single star, likely accompanied by some mass loss (Glebbeek et al. 2013). Since the system is no longer a binary, this does not lead to an X-ray binary phase nor a DNS system. Alternatively, the hydrogen-rich envelope is ejected before the two stars coalesce, leaving the stripped primary orbiting its companion (a main-sequence star in most cases). The result is thus similar to the case of stable mass transfer, except that orbital widening is avoided, and the orbit may have decayed significantly during the CE evolution (Ivanova et al. 2013). This is possibly the formation channel to short-orbit ($\lesssim 10\,\mathrm{days}$) HMXBs, including the RLO-HMXBs (Section 3.2). Also, as a result of the short duration of the CE inspiral, there is no expected spin-up of the main-sequence component. For further discussions of the CE phase, see Sections 1.5.1 and 4.2.

In the vast majority of cases, X-ray binary progenitors will consist of a star that is stripped of nearly all of its hydrogen-rich envelope and a main-sequence star companion. Binaries that avoid any mass transfer until the first SN are generally so wide that they would easily break up at this stage (Section 6), unless the stars evolve chemically homogeneously (e.g., de Mink et al. 2009; Marchant et al. 2016), in which case the orbit remains tight throughout the entire evolution.

Although depending on the initial parameters and the adopted physics of internal mixing, some models predict that the non-stripped component (i.e., the secondary star, the accretor from the first mass-transfer phase) will evolve to the SN stage first (Pols 1994; Wellstein et al. 2001). In such binaries—if subsequent reverse mass transfer could be avoided—the orbit would be more likely to be disrupted, since now the secondary star is more massive than its companion. In the more common situation of the stripped star evolving to the SN stage first (Figure 1), it is the less massive star exploding, and the binary would therefore not break up unless the NS receives a large kick at birth (Section 6).

We note that this picture, which is expected to apply to primary masses below about $40\,{M}_{\odot }$, may differ at higher masses. Above this threshold, the star radiates close to the Eddington limit, which may lead to loosely bound extended envelopes even in main-sequence stars (Sanyal et al. 2015) or to envelope ejection without the assistance of a companion as a luminous blue variable (Vanbeveren 1991) or yellow hypergiant star (Kourniotis et al. 2017). The consequences of these effects have not yet been well explored, but will also affect only a small fraction of the massive systems evolving to HMXBs.

For single stars at solar metallicity, the initial critical ZAMS mass to undergo a SN and produce a NS is $8\mbox{--}12\,{M}_{\odot }$ (Langer 2012 and references therein). This critical mass depends on the amount of convective core overshooting and the metallicity content of the star. For larger amounts of convective overshooting, and for lower metallicities, the threshold mass decreases significantly. In close binaries, this limit depends on the initial system parameters. Especially for mass donors in the initially closest systems, the mass limit at solar metallicity can be as high as $15\,{M}_{\odot }$ (Wellstein et al. 2001). On the other hand, as discussed by Podsiadlowski et al. (2004), electron capture SNe may occur in close binaries from stars with an initial mass below $8\,{M}_{\odot }$ because the second dredge-up episode, which reduces the helium core, can be avoided once the hydrogen envelope is lost.

HMXBs are composed of a compact object (NS or BH) orbiting a luminous and massive ($\gtrsim 10\,{M}_{\odot }$), early spectral type companion star. We can distinguish between three types of HMXBs according to their accretion process: (i) Be-star X-ray binaries (BeHMXBs), (ii) supergiant X-ray binaries (sgHMXBs), and (iii) Roche-lobe overflow systems (RLO-HMXBs). For a review of HMXBs, see, e.g., Chaty (2013) and Walter et al. (2015). We now discuss each of these three types in more detail.

3.2.1. BeHMXBs

BeHMXBs host a main-sequence donor star of spectral type B0–B2e III/IV/V, a rapid rotator surrounded by a circumstellar (so-called "decretion") disk of gas, as seen by the presence of a prominent Hα emission line. This disk is created by a low-velocity/high-density stellar wind of $\sim {10}^{-7}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (Lamers & Cassinelli 1999). Transient and bright (Type I) X-ray outbursts periodically occur each time the compact object (usually a NS in a wide and eccentric orbit) approaches periastron and accretes matter from the decretion disk (see Charles & Coe 2006 and references therein). These systems exhibit a correlation between NS spin and orbital period, as shown by their location from the lower center to the upper right of the Corbet diagram (see Figure 4), due to the efficient transfer of angular momentum at each periastron passage: rapidly spinning NSs correspond to short orbital period systems, and slowly spinning NSs to long orbit systems. The BeHMXBs are generally believed to be in spin equilibrium, meaning that the outer edge of the NS magnetosphere (defining the inner edge of the truncated accretion disk of the NS) rotates with Keplerian velocity (Davidson & Ostriker 1973; Stella et al. 1986; Waters & van Kerkwijk 1989). One can notice two outliers, however, SAX J2103.5+4545 and 1A1118-615, to this trend among BeHMXBs. Reig et al. (2004) and Staubert et al. (2011) explain these systems through long episodes of quiescence, presumably due to the lack of sufficient circumstellar matter in the decretion disk of the Be-star, whereby accretion mainly takes place from the stellar wind provided by the Be-star and thus causing the pulsars to spin down to the (longer) equilibrium period that is expected for wind-fed systems. Apart from MWC 656 possibly hosting a BH (Casares et al. 2014), the vast majority of (if not all other) known BeHMXBs seem to host a NS.

Zoom In Zoom Out Reset image size

3.2.2. sgHMXBs

3.2.3. RLO-HMXBs

A total of 114 HMXBs are reported in Liu et al. (2006), and 117 in Bird et al. (2016). By cross-correlating both catalogs, we compute that they share 79 HMXBs in common, although six of these common sources are actually now assigned to other types in the Bird et al. (2016) catalog. Adding new identifications by Coleiro & Chaty (2013) and Fortin et al. (2017), we find that the total number of HMXBs currently known in our Galaxy amounts to 167. These sources thus represent $\sim 43 \%$ of the total number of high-energy binary systems, i.e., adding all known LMXBs and HMXBs reported in Liu et al. (2007), Coleiro & Chaty (2013), Bird et al. (2016), and Fortin et al. (2017), respectively. Among the 167 HMXBs, there are 70 firmly identified BeHMXBs, 35 sgHMXBs (including 12 SFXTs), 3 RLO-HMXBs, and 59 HMXBs of unidentified nature (Fortin et al. 2017). Thus, in relative fractions, the subpopulations of HMXBs are 42% BeHMXBs, 21% sgHMXBs (7% SFXTs), 2% RLO-HMXBs, and 35% unidentified HMXBs.

Contrary to LMXBs (which have lost traces of their precise birth), HMXBs are concentrated in the Galactic plane, toward tangential directions of the Galactic arms, rich in star-forming regions and stellar complexes (Coleiro & Chaty 2013 and references therein). The distribution of HMXBs in any galaxy is thus a good indicator of star formation activity, and their collective X-ray luminosity has been used to compute the star formation rate of their host galaxy (Grimm et al. 2003). Coleiro & Chaty (2013) have correlated the position of Galactic HMXBs (including BeHMXBs and sgHMXBs) with the position of star-forming complexes (SFCs), and showed that HMXBs are clustered within $0.3\,\mathrm{kpc}$ of an SFC, with an intercluster distance of $1.7\,\mathrm{kpc}$, thus showing that HMXBs remain close to their birth place (see also Bodaghee et al. 2012, who found a slightly larger average offset between HMXBs and OB associations, but still consistent within the error bars of Coleiro & Chaty 2013). As we shall discuss later, in connection to the kinematics of DNS systems (Section 6), this is an important finding for constraining the magnitude of the kick velocity imparted on the resulting NS formed in the first SN explosion. In addition, Coleiro & Chaty (2013) showed that the HMXB distribution is offset by an age difference of $\sim {10}^{7}\,\mathrm{yr}$ with respect to the spiral arms, corresponding to the delay between the star birth and the HMXB formation. By taking into account the galactic arm rotation, they could derive parameters such as the age, migration distance, and systemic velocities for a number of BeHMXBs and sgHMXBs.

In Figure 4, we have plotted an "extended Corbet diagram," which includes the locations of the known DNS systems. An interesting question is how the different types of HMXBs and DNS systems are connected in an evolutionary context.

As the massive star in the different types of HMXBs evolves, there will be a point in most systems (except for those with very long orbital periods) when it starts to fill its Roche lobe and transfer mass by RLO. For NS systems, the ensuing mass transfer is generally unstable because of the large mass ratio which causes the orbit to shrink efficiently (e.g., Tauris & van den Heuvel 2006 and references therein). An additional effect that strongly contributes to unstable mass transfer is the Darwin instability¹⁶ (Darwin 1879). This leads to runaway accretion onto the NS and the build-up of a CE engulfing the whole system (see Sections 1.5.1 and 4.2). The further evolution of the system then depends on whether the CE can be ejected and on the fate of the NS inside the massive envelope. The following discussion assumes that the NS does not experience hypercritical accretion inside the CE, in which case it could potentially be converted to a BH (Chevalier 1993), which requires a different channel for forming DNSs (see the more detailed discussion in Section 4.2).

3.4.1. Close-orbit HMXBs and Formation of TŻOs

3.4.2. Wide-orbit HMXBs and the Formation of DNS Systems

The recent discovery that a number of ultra-luminous X-ray sources (ULXs: M82 X-2, NGC 7793 P13, and NGC 5907) are pulsating NSs (Bachetti et al. 2014; Fürst et al. 2016; Israel et al. 2017a, 2017b) has led to the suggestion that ULXs might be HMXBs (Ekşi et al. 2015; Kluźniak & Lasota 2015; Shao & Li 2015; Mushtukov et al. 2015), including BeHMXBs (Karino & Miller 2016; Christodoulou et al. 2017), and which could therefore suggest that some ULXs might potentially lead to the production of DNS systems.

The Chandra discovery of a population of ULXs with lifetimes less than about 10 Myr in the Cartwheel galaxy suggests indeed that ULXs could be related to HMXBs (King 2004). However, these ULX sources were not detected to be pulsating and might therefore host accreting BHs. If they host BHs (which are more massive than NSs), their orbital evolution will be much more stable.

ULXs have also been suggested to be intermediate-mass X-ray binaries (IMXBs; King & Lasota 2016; Karino 2016; Israel et al. 2017a), which have typical donor stars of $3\mbox{--}8\,{M}_{\odot }$. Theoretically, it has been established that the orbital evolution of IMXBs is dynamically stable for a range of initial donor star masses and orbital periods (Tauris et al. 2000; Podsiadlowski et al. 2002; Shao & Li 2012), and that these systems therefore avoid the formation of a CE (King & Begelman 1999). During the X-ray phase, IMXBs might mimic systems like SS 433 (King et al. 2000) and eventually they may evolve into systems like Cyg X-2 (King & Ritter 1999; Podsiadlowski & Rappaport 2000). The final product of such IMXBs, however, is not a DNS system but instead a recycled NS orbited by a (typically CO) WD companion (Tauris et al. 2000, 2011; Podsiadlowski et al. 2002; Lin et al. 2011).

To investigate whether ULXs are indeed IMXBs or HMXBs (and in the latter case, potential progenitors of DNS systems), we have analyzed the orbital evolution of such systems. In Figure 5, we have plotted the relative change in orbital separation ($a/{a}_{0}$) as a function of the decreasing donor star mass, M₂ for a HMXB (initially ${M}_{2}=12\,{M}_{\odot }$) and an IMXB (initially ${M}_{2}=6\,{M}_{\odot }$). To calculate these tracks of non-conservative RLO, we applied the so-called isotropic re-emission model (Bhattacharya & van den Heuvel 1991; van den Heuvel 1994). By integrating the orbital angular momentum balance equation, one can derive the analytical expression (Tauris 1996)

$$\begin{eqnarray}&&\displaystyle \frac{a}{{a}_{0}}={\left(\displaystyle \frac{{q}_{0}(1-\beta )+1}{q(1-\beta )+1}\right)}^{\tfrac{3\beta -5}{1-\beta }}\left(\displaystyle \frac{{q}_{0}+1}{q+1}\right){\left(\displaystyle \frac{{q}_{0}}{q}\right)}^{2},\end{eqnarray} \tag{ 2 }$$

where a₀ and a refer to the orbital separation before and during RLO, respectively, and where q₀ and q represent the mass ratios at these two epochs. The parameter β represents the (assumed constant) fraction of transferred material that is lost from the vicinity of the NS. Here we assume that the wind mass loss rate is negligible compared to the rate at which material is transferred via L1 during RLO, and we disregard the possibility of formation of a circumbinary disk.

Zoom In Zoom Out Reset image size

Given that the thermal timescale mass-transfer rate resulting from RLO in HMXB and IMXB systems is much higher than the Eddington limit for an accreting NS ($\sim 1.8\times {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$), we expect evolution close to $\beta \approx 1$. There is evidence from strong spin-up torques that NSs in ULXs are able to accrete substantially above the Eddington limit (even on the order of $\sim 100\,{\dot{M}}_{\mathrm{Edd}}$; Israel et al. 2017a). However, given that IMXBs and HMXBs immediately reach values of $| {\dot{M}}_{2}| \gt {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ after the onset of RLO, the estimate of $\beta \gt 0.9$ will likely still remain valid.

From Figure 5, it can be seen that the orbital shrinkage is substantially larger for a HMXB system compared to an IMXB system by more than an order of magnitude. Therefore, as discussed previously, HMXBs are bound to become dynamically unstable upon efficient accretion (delayed dynamical instability (DDI); Hjellming & Webbink 1987), unless some (most likely unrealistic) very large value of the wind mass loss rate exceeding the rate at which matter is transferred via L1 during RLO (i.e., $\alpha \gt 0.9$, Fragos et al. 2015), or other unconventional types of applied angular momentum losses (however, see also Pavlovskii et al. 2017), is assumed.

The (thermal) timescale on which such a HMXB becomes dynamically unstable after initiating RLO (and presumably forming a CE) is of the order ${10}^{2}\mbox{--}{10}^{4}\,\mathrm{yr}$ (Savonije 1978), depending on the orbital period at RLO. In comparison, the IMXB systems are typically able to provide mass-transfer rates $| {\dot{M}}_{2}| \gt {10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ for a few ${10}^{5}\mbox{--}{10}^{6}\,\mathrm{yr}$ (e.g., Tauris et al. 2000, 2011), followed by RLO at a lower rate for some ${10}^{6}\mbox{--}{10}^{7}\,\mathrm{yr}$, depending on the orbital period. Therefore, based on the above arguments and timescale estimates, we find it more plausible that any detected NS ULX system is an IMXB source compared to a HMXB source (although the latter case cannot be ruled out). Hence, we find that the recently observed pulsating NS ULXs are unlikely to be progenitors of DNS systems. The main points of the above discussions are shown in Table 4.

Table 4.  Linking IMXBs, HMXBs, NS ULXs, and DNS Systems

${M}_{2}/{M}_{\odot }$	Stability	NS ULX source	Final outcome
$\lesssim 5$	stable RLO	yes	NS + He/CO WD
6–9	$\mathrm{DDI}\to \mathrm{CE}$	?	NS + CO/ONeMg WD
$\gtrsim 10$	CE	no	DNS

Note. The division between IMXBs and HMXBs is somewhat unclear but typically at a donor star mass of ${M}_{2}\simeq 8\mbox{--}10\,{M}_{\odot }$. A DNS system is only produced from a HMXB evolving through a CE phase if the initial orbital period is wide enough (see Section 3.4.2).

Download table as:  ASCIITypeset image

We have recently been notified that optical monitoring of the NS ULX system NGC 7793 P13 (Motch et al. 2014) has revealed a $\sim 20\,{M}_{\odot }$ donor star and an orbital period of ∼64 days. If this interpretation is correct, it is a very puzzling result and it challenges current understanding of long-term stability of such HMXBs (unless this source initiated RLO very recently and was fortunate to be discovered in the short time span before it becomes dynamically unstable). If the system in the future evolves through a CE phase, then we expect it to coalesce and produce a TŻO, as explained in Section 3.4).

To summarize our analysis of DNS progenitor systems in this section, we conclude that only wide-orbit HMXBs (possibly only a few, if any, of the observed Galactic HMXBs) will eventually survive the CE evolution and potentially produce DNS systems. NS ULX systems are possibly IMXBs and, if so, not massive enough to be progenitors of DNS systems.

The stellar wind mass-loss rate of early type stars is subject to large uncertainties. The observational and theoretical estimated values depend on the modeling of the terminal wind velocity, clumping effects, metallicity, and absorption properties in the local environment (Kudritzki & Puls 2000; Vink et al. 2001; Repolust et al. 2004; Crowther et al. 2006; Mokiem et al. 2007; Markova & Puls 2008; Langer 2012; Šurlan et al. 2013; Chaty 2013; Falanga et al. 2015). In the following, we will begin by assuming a wind mass-loss rate of a typical OB-star companion in a HMXB of $| {\dot{M}}_{2,\mathrm{wind}}| \simeq {10}^{-6}\mbox{--}{10}^{-5}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (assuming that the effect of X-ray irradiation on the wind mass-loss rate of the companion star is negligible). The wind velocity can be approximated by the escape velocity at the stellar surface, ${v}_{\mathrm{wind}}\simeq {v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{2}/{R}_{2}}$, where M₂ and R₂ are the companion star mass and radius, and G is the constant of gravity. Moreover, we assume the Bondi–Hoyle-like accretion to be valid, given the highly supersonic flow of the stellar wind from the companion star. The gravitational capture radius of the accreting NS is then given by ${R}_{\mathrm{acc}}=2{{GM}}_{\mathrm{NS}}/({v}_{\mathrm{rel}}^{2}+{c}_{s}^{2})$, where ${c}_{s}=\sqrt{\gamma P/\rho }\approx 11\,\sqrt{T/({10}^{4}\,K)}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ is the local sound speed of the ambient medium (γ is the adiabatic index, P is the pressure, ρ is the mass density, and T is the temperature). The relative velocity between the wind and the NS (${{\boldsymbol{v}}}_{\mathrm{rel}}={{\boldsymbol{v}}}_{\mathrm{wind}}+{{\boldsymbol{v}}}_{\mathrm{orb}}$) is usually dominated by the wind velocity compared to the orbital velocity, i.e., ${v}_{\mathrm{wind}}\gg {v}_{\mathrm{orb}}$ (except for very tight and/or highly eccentric binaries). Similarly, we have ${v}_{\mathrm{wind}}\gg {c}_{s}$. The resulting accretion rate onto the NS, ${\dot{M}}_{\mathrm{NS}}\approx \pi {R}_{\mathrm{acc}}^{2}\rho \,{v}_{\mathrm{wind}}$, can be combined with the continuity equation $\rho =| {\dot{M}}_{2,\mathrm{wind}}| /(4\pi {a}^{2}\,{v}_{\mathrm{wind}})$ to yield

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{NS}}\approx \displaystyle \frac{{({{GM}}_{\mathrm{NS}})}^{2}\,| {\dot{M}}_{2,\mathrm{wind}}| }{{a}^{2}{v}_{\mathrm{wind}}^{4}},\end{eqnarray} \tag{ 3 }$$

where a is the orbital separation. An estimate from this equation indicates that, in wind-driven scenarios, ${\dot{M}}_{\mathrm{NS}}$ is typically ${10}^{3}-{10}^{4}$ times lower than $| {\dot{M}}_{2,\mathrm{wind}}|$ (for helium star winds, the ratio is even lower by an order of magnitude; see Section 4.3). Hence, we adopt a typical NS accretion rate of ${\dot{M}}_{\mathrm{NS}}\approx {10}^{-9}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, corresponding to an X-ray luminosity of the order ${L}_{{\rm{x}}}\approx {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$.

To estimate the amount of material accreted by the NS during wind accretion in the HMXB stage, we then simply multiply this rate with the expected lifetime of the HMXB companion star during which NS wind accretion is active: ${\rm{\Delta }}{M}_{\mathrm{NS}}={\dot{M}}_{\mathrm{NS}}\cdot {\rm{\Delta }}{t}_{2}$. However, an analysis of the X-ray luminosity function of Galactic HMXBs (Grimm et al. 2002) shows that only ∼10% of the sources have ${L}_{{\rm{x}}}\geqslant {10}^{37}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$. This fact probably reflects a limited amount of time in which a HMXB companion star is able to deliver a sufficiently high wind mass-loss rate and/or our initially assumed (average) wind mass-loss rate, $| {\dot{M}}_{2,\mathrm{wind}}|$, is too large. During most of the (main-sequence) lifetime of the HMXB donor star, the value of $| {\dot{M}}_{2,\mathrm{wind}}|$ is often lower than ${10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. Indeed, from a study of 10 eclipsing HMXBs, Falanga et al. (2015) list typical values of $| {\dot{M}}_{2,\mathrm{wind}}| \approx {10}^{-7}\mbox{--}{10}^{-6}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, based on observations of their average NS accretion rates. Hence, by adopting a value of ${\rm{\Delta }}{t}_{2}\sim {10}^{7}\,\mathrm{yr}$ and slightly smaller values of $| {\dot{M}}_{2,\mathrm{wind}}|$ and ${\dot{M}}_{\mathrm{NS}}$, we finally obtain ${\rm{\Delta }}{M}_{\mathrm{NS}}\approx {\rm{a}}\,\mathrm{few}\,\times {10}^{-3}\,{M}_{\odot }$ for wind accretion during the HMXB phase.

The structured nature of the stellar wind of a supergiant star (i.e., a wind where cool dense clumps are embedded in a photoionized gas) makes it difficult to derive an estimate of its mass-loss rate that is more accurate than a factor of a few (e.g., Martínez-Núñez et al. 2017). Another caveat is that the presence of a NS causes hydrodynamical interactions of the wind flow in close binaries (Manousakis et al. 2012). However, it is also possible that a strong reduction of the long-term mass accretion rate is produced by the effect of a strong magnetic field and/or a fast spin of the NS, which could lead to magnetic and/or centrifugal barriers, respectively (Bozzo et al. 2008, 2016). Such complications are not taken into account in our simplified Bondi–Hoyle-like treatment of the accretion process (see also Section 4.6).

The question of accretion onto a NS undergoing spiral-in during a CE event has been the subject of long debate. From a theoretical point of view, it has been suggested that a NS in a CE might suffer from hypercritical accretion (Chevalier 1993) leading to its collapse into a BH. Thus, to explain the observed DNS systems, an alternative scenario without CE evolution was invented: the double core scenario (Brown 1995; Dewi et al. 2006). In this scenario, two stars with an initial mass ratio close to unity evolve in parallel and reach the giant stage roughly at the same time. Therefore, when the CE forms, it will embed both stars in their giant stages (or as a giant star and a helium core), thereby avoiding the formation of a CE with a NS. This scenario, however, only works for the evolution of two stars with a mass ratio close to 1 and thus cannot explain the observations of tight binaries with, for example, a NS orbited by a WD companion, e.g., PSR B0655+64 (Damashek et al. 1982; Kulkarni 1986), PSR J1802−2124 (Ferdman et al. 2010), and PSR J1952+2630 (Lazarus et al. 2014). Such systems show clear evidence of having evolved through a stage with a NS embedded in a CE (van den Heuvel & Taam 1984; Tauris et al. 2012).

Recent theoretical works by MacLeod & Ramirez-Ruiz (2015a, 2015b) seem to resolve this paradox. From their hydrodynamical simulations, they find that a NS embedded in a CE only accretes a very modest amount of material during its spiral-in as a result of a density gradient across its accretion radius, which strongly limits accretion by imposing a net angular momentum on the flow around the NS. MacLeod & Ramirez-Ruiz (2015b) find an upper limit of ${\rm{\Delta }}{M}_{\mathrm{NS}}\lt 0.1\,{M}_{\odot }$. However, they remark in their paper that they make several approximations that likely lead their calculations of accreted mass in a CE to be an overestimate, including an assumption of a static structure for the CE, effective hypercritical accretion, and cooling by neutrinos, as well as propagation of accreted mass all the way down to the shock radius.

Observational support for such a small amount of NS accretion during the CE phase is given by the four DNS systems (Table 2) in which the first-born NS (the recycled component) has a mass of less than $1.38\,{M}_{\odot }$, and also by PSR J1802−2124, where the NS mass is estimated to be $\sim 1.24\pm 0.11\,{M}_{\odot }$ (Ferdman et al. 2010). If all of these NSs had accreted of the order $0.1\,{M}_{\odot }$, they would need to have been born with a mass of ${M}_{\mathrm{NS}}\lt 1.28\,{M}_{\odot }$. While this is still possible, it would be unexpected for the mass of the first-born NS in so many systems. Furthermore, the NS is also expected to accrete material in the subsequent phases discussed below, which would require even smaller birth masses in these cases. In the following, we take ${\rm{\Delta }}{M}_{\mathrm{NS}}=0.01\,{M}_{\odot }$ to be a reasonable estimate for the upper limit of the amount of mass accreted by a NS during the CE phase, although we cannot fully rule out the possibility of further accretion.

The recipe for calculating the amount of mass accreted by the NS from the fast wind of its naked helium (or Wolf–Rayet, WR) star companion is similar to that given in Section 4.1, i.e., ${\rm{\Delta }}{M}_{\mathrm{NS}}\simeq | {\dot{M}}_{\mathrm{He},\mathrm{wind}}| \cdot {f}_{\mathrm{acc}}\cdot {\rm{\Delta }}{t}_{\mathrm{He}}$, where ${\rm{\Delta }}{t}_{\mathrm{He}}$ is the lifetime of the helium star, $| {\dot{M}}_{\mathrm{He},\mathrm{wind}}|$ is its mass-loss rate, and ${f}_{\mathrm{acc}}\,={\dot{M}}_{\mathrm{NS}}/| {\dot{M}}_{\mathrm{He},\mathrm{wind}}|$ is the wind accretion efficiency. It should be noted that helium stars with mass $\lesssim 8\,{M}_{\odot }$ are never seen as WR stars, meaning that strong WR winds occur only in helium stars $\gtrsim 8\,{M}_{\odot }$ (Crowther 2007). Such massive helium stars, however, rarely end up as progenitors of DNS systems. From the stellar evolution modeling of NS–helium star binaries (Tauris et al. 2015), we find that typically ${f}_{\mathrm{acc}}\simeq {\rm{a}}\,\mathrm{few}\,\times \,{10}^{-4}$ (see also Israel 1996), and we obtain values of ${\rm{\Delta }}{M}_{\mathrm{NS}}\lt 4\times {10}^{-4}\,{M}_{\odot }$ when integrating throughout the wind accretion phase of the NS–helium star binaries.

The amount of material that can be accreted by the NS as a result of Case BB RLO is limited by the short duration of this mass-transfer phase ($\lt {10}^{5}\,\mathrm{yr}$) combined with an Eddington-limited accretion rate onto the NS, which is roughly given by (cgs units in the central term)

$$\begin{eqnarray}&&{\dot{M}}_{\mathrm{Edd}}=\displaystyle \frac{4\pi c\,{R}_{\mathrm{NS}}}{0.2(1+X)}\,\simeq \,3.0\times {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}\end{eqnarray} \tag{ 4 }$$

for the accretion of helium (X = 0). Applying detailed binary stellar evolution modeling, Tauris et al. (2015) systematically investigated a wide parameter space of initial orbital periods and helium star masses to calculate the amount of mass accreted by the NS. For binaries leading to the formation of DNS systems, they found ${\rm{\Delta }}{M}_{\mathrm{NS}}=5\times {10}^{-5}\mbox{--}3\times {10}^{-3}\,{M}_{\odot }$, where ${\rm{\Delta }}{M}_{\mathrm{NS}}$ is systematically decreasing with increasing values of ${P}_{\mathrm{orb}}$. These values are in agreement with previous studies of Case BB RLO (Dewi et al. 2002; Dewi & Pols 2003; Ivanova et al. 2003; Tauris et al. 2012; Lazarus et al. 2014).

Besides increasing the mass of the accreting NS star, Case BB RLO also serves to transfer angular momentum and thus spin up this first-born NS to become observable as a (mildly) recycled radio pulsar. Based on standard assumptions of the torque acting on the NS magnetosphere at the inner edge of the accretion disk, combined with observational data of recycled pulsars in DNS systems, Tauris et al. (2012) and Lazarus et al. (2014) concluded that Case BB RLO should allow for NS accretion rates a factor of at least 2–3 above the standard value of ${\dot{M}}_{\mathrm{Edd}}$ (see above). This result seems to be necessary to explain the fastest spins of recycled pulsars with NS or massive WD companions (which are also believed to have evolved via post-CE Case BB RLO). An X-ray binary can circumvent the value of ${\dot{M}}_{\mathrm{Edd}}$, for example, due to a strong B-field of the NS (see discussions in Section 5.2). Hence, we conclude that Case BB RLO may, in a few cases, result in accretion of up to $(6\mbox{--}9)\times {10}^{-3}\,{M}_{\odot }$.

In case the secondary star explodes at a very close distance to the first-born NS, the latter may, in principle, accrete material from the SN ejecta (Fryer et al. 2014). However, in order for this to happen, the evolution leading to the second SN event has to be extremely fine-tuned with respect to both the preceding CE phase and the evolutionary stage of the secondary star. The system would have to almost, but not quite, merge during the CE and at the same time ensure that the secondary star explodes before the system merges as a result of GW radiation. Although one cannot rule out that such a fine-tuning might be possible in a few extremely rare cases—leading to orbital periods of 2–10 minutes at the moment of explosion—such a scenario can safely be ignored for DNS progenitors in general and thus ${\rm{\Delta }}{M}_{\mathrm{NS}}\ll {10}^{-3}\,{M}_{\odot }$.

Although the HMXB donor star (the progenitor of the second-born NS) provides the material for potential accretion onto the first-born NS in all of the above-mentioned phases, it is not certain how much of that material transferred will actually be accreted. For example, some of it may be lost from the system due to ejector or propeller effects, accretion disk instabilities, and direct irradiation of the donor's atmosphere from the pulsar (e.g., Illarionov & Sunyaev 1975; van Paradijs 1996; Dubus et al. 1999; Romanova et al. 2009). Evidence of low NS accretion efficiencies in LMXB binaries is seen in several binary pulsar systems with fully recycled MSPs but with relatively small NS masses (Tauris & Savonije 1999; Antoniadis et al. 2012, 2017). An extreme example is PSR J1918−0642, which is an MSP with $P=7.6\,\mathrm{ms}$, in orbit with a low-mass He WD and ${P}_{\mathrm{orb}}=10.9\,\mathrm{days}$, and yet it has a mass of only ${1.18}_{-0.09}^{+0.10}\,{M}_{\odot }$ (Fonseca et al. 2016).

As discussed in Section 3.4, only wide-orbit HMXBs are expected to survive CE evolution and produce DNS systems (e.g., Taam 1996). For such HMXBs, both ${\rm{\Delta }}{M}_{\mathrm{NS}}$ and the efficiency of angular momentum transfer are even more uncertain. Moreover, there is evidence of quite a variety in NS accretion efficiencies from the location of the different populations of HMXBs in the Corbet diagram (Waters & van Kerkwijk 1989; Chaty 2013; Li et al. 2016). It is therefore likely that the estimated values of ${\rm{\Delta }}{M}_{\mathrm{NS}}$ at various stages discussed above may represent overestimates due to magnetospheric and accretion disk effects, which were not taken into account.

Adding the expected amounts of material accreted by the first-born NS from each of the above-mentioned stages leads to a total accumulated amount of ${\rm{\Delta }}{M}_{\mathrm{NS}}\lt 0.02\,{M}_{\odot }$. The exact amount depends on CE physics and magnetospheric conditions, i.e., the B-field and the spin evolution of the NS. The spin-up, on the other hand, is expected to be most efficient during Case BB RLO where an accretion disk is always present (unlike the situation in a CE; Murguia-Berthier et al. 2017), leading to a strong and stable accretion torque acting on the NS. Hence, the spin period of the recycled pulsar is primarily determined by the amount of accreted material during Case BB RLO.

The minimum amount of mass accreted to recycle a pulsar to reach a certain equilibrium spin period is estimated to be roughly given by (see discussions in Sections 4.3 and 4.4 in Tauris et al. 2012)¹⁷

$$\begin{eqnarray}&&{\rm{\Delta }}{M}_{\mathrm{eq}}\simeq 0.22\,{M}_{\odot }\,\displaystyle \frac{{({M}_{\mathrm{NS}}/{M}_{\odot })}^{1/3}}{{P}_{\mathrm{ms}}^{4/3}}.\end{eqnarray} \tag{ 5 }$$

From the fact that the observed spin periods of the (mildly) recycled NSs in DNS systems are in the interval $23\mbox{--}185\,\mathrm{ms}$ (see Table 2), we notice that only a very modest amount, ${\rm{\Delta }}{M}_{\mathrm{NS}}=(0.2\mbox{--}4)\times {10}^{-3}\,{M}_{\odot }$, is in principle needed to explain these spin periods. Interestingly enough, these values of ${\rm{\Delta }}{M}_{\mathrm{NS}}$ are in fine agreement with the amount of accreted material derived above.

According to the mass-transfer modeling of Tauris et al. (2015), the fastest spinning recycled DNS pulsars (outside of globular clusters) are expected to possess spin periods down to $\sim 11\,\mathrm{ms}$. Hence, no DNS systems are expected to be discovered with a spin period below this value. This argument strengthens the hypothesis that the globular cluster DNS J1807−2500B, which has a spin period of only 4.2 ms, is the outcome of a dynamical encounter event (Lynch et al. 2012). This NS was most likely spun-up in a LMXB system and where the low-mass companion (WD) was subsequently exchanged with a more massive ($1.21\,{M}_{\odot }$) compact object to produce the presently observed DNS system.

A caveat in the above discussion is the possibility of long-term accretion at a rate ${\dot{M}}_{\mathrm{NS}}\gg {\dot{M}}_{\mathrm{Edd}}$. Such a high accretion rate is apparently evident in a ULX system with a pulsating NS accretor (Israel et al. 2017a). If such a ULX system is indeed a HMXB (rather than an IMXB, as we argue in favor of in Section 3.5), and if such a system could avoid a delayed dynamical instability and actually form a DNS system, then it might be possible to produce a DNS system in the Galactic disk with a NS being fully recycled to a spin period of $1\mbox{--}2\,\mathrm{ms}$ (or even a sub-ms pulsar). A remaining puzzle, however, is how the strong B-field (if this is the reason behind a highly super-Eddington accreting NS) would remain strong without decaying (Bhattacharya 2002), given the vast amount of material accreted by such a NS. If a hypothetical post-ULX DNS system hosted a recycled pulsar with a large B-field, it would only be detectable for a short time due to the rapid loss of rotational energy. However, the second-born NS would be detectable as a normal pulsar.

The above analysis demonstrates the very limited accretion ($\lt 0.02\,{M}_{\odot }$) by the first-born NS during/after the HMXB phase, and we can hereby conclude that the distribution of NS masses in DNS systems (Figure 3), either the first-born (mildly recycled) NS or the second-born (young) NS, to a large degree reflects the birth distribution of NS masses in such observed systems. Two interesting questions are therefore: (i) whether the mass difference (roughly $\sim 0.1\,{M}_{\odot }$) between the first- and the second-born NSs can be understood from a stellar evolution (and/or a SN explosion) point of view, and (ii) how the birth masses of the NSs in DNS systems compare to the measured masses of NSs observed in X-ray binaries or in binary radio pulsar systems with a WD companion.

To answer the first question, it is of interest to consider studies of the pre-SN stripping effects of the exploding stars in close orbits (SNe Ib/c). Due to the presence of a non-degenerate companion star in the first SN, it is clear that the progenitor of the first-born NS (Yoon et al. 2010) is usually stripped significantly less than the progenitor of the second-born NS (Tauris et al. 2015). The reason is that in the second SN, the progenitor of the exploding star is stripped significantly deeper when the companion is a NS. This might give a natural explanation for the potential differences in the resulting NS masses. Note that we must be aware of potential selection effects. If a small NS mass is connected to a small kick (see Section 6.6), the systems with the smallest masses of the second-born NSs have a larger probability of surviving the second SN. This provides a selection effect favoring low secondary NS masses to be overrepresented among the DNS systems. In Section 6.9, we demonstrate, however, that this selection effect is rather limited in practice.

For the second question, we see in Figure 3 that while the lower end of the NS mass distribution matches fairly well between NSs in DNS and NS+WD systems, so far no massive NS ($\geqslant 1.70\,{M}_{\odot }$) has been found in a DNS system. However, we are dealing with relatively small number statistics and we know of NS masses in HMXBs (the progenitors of DNS systems) that may be as high as $\sim 2.0\,{M}_{\odot }$, e.g., Vela X-1 (Barziv et al. 2001; Falanga et al. 2015). Given that wide-orbit HMXBs are progenitors of DNS systems (Section 3.4), one would indeed expect the NS masses in these HMXBs to resemble the masses of the first-born NSs in DNS systems based on our presented evidence for very limited accretion during/after the HMXB phase. Unfortunately, the masses of NSs in HMXBs (Falanga et al. 2015) are not measured to the same precision as radio pulsar masses, which limits the possibility of a useful direct comparison.

In Figure 6, we have plotted the spin period of the recycled pulsars as a function of ${P}_{\mathrm{orb}}$ for all observed Galactic disk DNS systems. A fit to the raw data results in the following correlation (gray line):

$$\begin{eqnarray}&&P=44\,\mathrm{ms}\,{({P}_{\mathrm{orb}}/\mathrm{days})}^{0.26},\end{eqnarray} \tag{ 6 }$$

based on a linear regression in the ($\mathrm{log}{P}_{\mathrm{orb}},\mathrm{log}P$) plane and with a regression coefficient of ${R}^{2}=0.73$. However, as we shall demonstrate below, differences in the progenitor stars and their SN explosions result in a scatter in this diagram. Furthermore, the observed data of DNS systems do not necessarily reflect their birth properties of ${P}_{\mathrm{orb}}$ and P, since pulsars lose rotational energy with time, and for tight DNS systems, the orbits decay due to GW radiation (i.e., P increases and ${P}_{\mathrm{orb}}$ decreases; see Section 7). By taking these effects into account in a qualitative manner and placing particular weight on the very wide-orbit DNS 1930−1852 (Swiggum et al. 2015), we obtain the following simple empirical birth correlation (green line):

$$\begin{eqnarray}&&P\approx 36\pm 14\,\mathrm{ms}\,{({P}_{\mathrm{orb}}/\mathrm{days})}^{0.40\pm 0.10}.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

A corresponding theoretical $({P}_{\mathrm{orb}},P)$ correlation can also be derived for DNS systems and which can easily be understood qualitatively from binary stellar evolution arguments (Tauris et al. 2015). The wider the initial orbit of the DNS progenitor system (i.e., the NS–helium star system following the CE phase), the more evolved is the helium star by the time it fills its Roche lobe and initiates mass transfer. Hence, in a wide system, the helium star has less time to transfer material toward the accreting NS before it terminates its nuclear evolution and explodes in a SN. As a result, the wider the progenitor system is, the less efficient is the recycling of the first-born NS and the larger (slower) is its rejuvenated spin period. For example, helium star donors in wide-orbit systems may already be undergoing carbon-shell burning at the onset of Case BB RLO (strictly speaking Case BC RLO), whereby their remaining lifetime is less than a few thousand years before they collapse. Such a case is demonstrated in Figure 7. In the widest DNS systems, one would therefore expect to find NSs with $P\gt 100\,\mathrm{ms}$ and relatively large B-fields, i.e., marginally recycled pulsars (Tauris et al. 2015). The recently discovered DNS system PSR J1930−1852 (Swiggum et al. 2015) is an excellent example of such a case ($P=185\,\mathrm{ms}$, $\dot{P}=1.8\,\times {10}^{-17}\,{\rm{s}}\,{{\rm{s}}}^{-1}$, ${P}_{\mathrm{orb}}=45\,\mathrm{days}$, and eccentricity e = 0.40; see Figure 2). Its discovery is important for testing formation scenarios of such wide-orbit marginally recycled pulsars.

Zoom In Zoom Out Reset image size

Based on a few assumptions, one can derive a theoretical $({P}_{\mathrm{orb}},P)$ correlation as follows. From the detailed grid of Case BB RLO models by Tauris et al. (2015), we obtain the amount of material accreted by the NS, ${\rm{\Delta }}{M}_{\mathrm{NS}}$. Assuming the NS reaches equilibrium spin during Case BB RLO, we can use Equation (5) to estimate the spin period of the recycled pulsar, P, at the end of the mass-transfer phase (see also Dewi et al. 2005). The orbital period at the end of Case BB RLO, just prior to the (ultra-stripped) SN explosion, is also provided by the models in Tauris et al. (2015). Assuming, as a first-order approximation, that this second SN is symmetric (without any kick imparted on the newborn NS, i.e., w = 0; see Section 6 for discussions), we can use simple analytic prescriptions to calculate the post-SN orbital period, ${P}_{\mathrm{orb}}$, and the eccentricity (see Section 6.3).

Figure 8 shows examples of such theoretical $({P}_{\mathrm{orb}},P)$ correlations based on purely symmetric SNe. In this plot, we assumed in all cases an exploding star, which is a stripped helium star with an initial mass of $3.0\,{M}_{\odot }$. The parameters of the applied models from Tauris et al. (2015) are shown in Table 5. Three different blue curves are shown in Figure 8, depending on the assumed value for the Eddington accretion efficiency parameter, ${X}_{\mathrm{Edd}}$. This parameter is a measure of the maximum NS accretion rate in units of the classical Eddington limit, ${X}_{\mathrm{Edd}}\equiv {\dot{M}}_{\mathrm{NS}}^{\max }/{\dot{M}}_{\mathrm{Edd}}$. The larger the value of ${X}_{\mathrm{Edd}}$, the more efficient the recycling (more material is accreted by the NS), and the smaller the final spin period of the recycled pulsar. The quoted values of ${\rm{\Delta }}{M}_{\mathrm{NS}}$ in Tables 5 and 6 are based on an assumed accretion efficiency value of ${X}_{\mathrm{Edd}}=1$. However, as discussed in Section 4.4, it is possible that the amount of material accreted by the NS, ${\rm{\Delta }}{M}_{\mathrm{NS}}$, is larger by a factor of three, based on evidence for high accretion efficiency in PSR J1952+2630 (Lazarus et al. 2014).

Zoom In Zoom Out Reset image size

Table 5.  Properties of Case BB RLO Systems with ${M}_{\mathrm{He},{\rm{i}}}=3.0\,{M}_{\odot }$

${M}_{\mathrm{He},{\rm{i}}}$	${P}_{\mathrm{orb},{\rm{i}}}$	${M}_{\mathrm{He},{\rm{f}}}$	${M}_{\mathrm{core},{\rm{f}}}$	${P}_{\mathrm{orb},{\rm{f}}}$	${\rm{\Delta }}{M}_{\mathrm{NS}}$	${M}_{\mathrm{NS},2}$
(M⊙)	(days)	(M⊙)	(M⊙)	(days)	(M⊙)	(M⊙)
3.0	0.08	1.49	1.45	0.052	2.3e−3	1.31
3.0	0.10	1.58	1.50	0.063	2.0e−3	1.35
3.0	0.50	1.73	1.56	0.311	7.2e−4	1.40
3.0	2.00	1.80	1.57	1.23	4.8e−4	1.40
3.0	5.00	1.86	1.57	4.61	4.1e−4	1.40
3.0	10.0	2.00	1.57	9.51	2.3e−4	1.40
3.0	20.0	2.17	1.58	20.0	1.4e−4	1.41
3.0	50.0	2.32	1.58	53.2	6.3e−5	1.41
3.0	80.0	2.37	1.57	87.1	5.2e−5	1.40
3.0	100.0	2.38	1.57	109.7	4.8e−5	1.40

Note. Data taken from the calculations of Tauris et al. (2015) based on an initial NS accretor of mass ${M}_{\mathrm{NS}}=1.35\,{M}_{\odot }$ and initial orbital periods in the interval ${P}_{\mathrm{orb},{\rm{i}}}=0.08\mbox{--}100\,\mathrm{days}$. The quantities in the columns (left to right) are the initial helium star mass (${M}_{\mathrm{He},{\rm{i}}}$) and orbital period (${P}_{\mathrm{orb},{\rm{i}}}$), the final mass (and metal core mass) of the exploding helium star (${M}_{\mathrm{He},{\rm{f}}}$ and ${M}_{\mathrm{core},{\rm{f}}}$), the final orbital period prior to the SN explosion (${P}_{\mathrm{orb},{\rm{f}}}$), the amount of mass accreted by the NS (${\rm{\Delta }}{M}_{\mathrm{NS}}$) assuming ${X}_{\mathrm{Edd}}=1.0$, and the estimated mass of the newborn secondary NS (${M}_{\mathrm{NS},2}$)—see Section 5.2.

Download table as:  ASCIITypeset image

Table 6.  Properties of Various Case BB RLO Systems

${M}_{\mathrm{He},{\rm{i}}}$	${P}_{\mathrm{orb},{\rm{i}}}$	${M}_{\mathrm{He},{\rm{f}}}$	${M}_{\mathrm{core},{\rm{f}}}$	${P}_{\mathrm{orb},{\rm{f}}}$	${\rm{\Delta }}{M}_{\mathrm{NS}}$	${M}_{\mathrm{NS},2}$
(M⊙)	(days)	(M⊙)	(M⊙)	(days)	(M⊙)	(M⊙)
2.8	0.10	1.43	1.39	0.079	2.7e−3	1.26
3.0	0.10	1.58	1.50	0.063	2.0e−3	1.35
3.5	0.10	1.88	1.76	0.048	1.2e−3	1.56
4.0	0.10	2.12	1.96	0.048	9.4e−4	1.71
6.0	0.10	2.37	2.15	0.064	6.7e−4	1.86
2.7	0.50	1.48	1.41	0.415	8.4e−4	1.27
3.0	0.50	1.73	1.56	0.311	7.2e−4	1.40
3.5	0.50	2.07	1.80	0.365	5.0e−4	1.59
2.6	2.00	1.46	1.37	1.78	8.3e−4	1.24
3.0	2.00	1.80	1.57	1.23	4.8e−4	1.40
3.5	2.00	2.39	1.81	1.78	2.6e−4	1.60
4.0	2.00	2.95	2.02	2.45	1.4e−4	1.76
2.6	20.0	1.56	1.37	19.3	6.3e−4	1.24
3.0	20.0	2.17	1.58	20.0	1.4e−4	1.41
3.2	20.0	2.67	1.70	23.8	2.8e−5	1.51
3.0	50.0	2.32	1.58	53.2	6.3e−5	1.41
3.0	80.0	2.37	1.57	87.1	5.2e−5	1.40

Note. Data from Tauris et al. (2015). See notes for Table 5.

Download table as:  ASCIITypeset image

Super-Eddington accretion by a factor of a few is easily obtained by a thick accretion disk with a central funnel (Abramowicz et al. 1980) or a reduction in the electron-scattering cross-section caused by a strong B-field (Basko & Sunyaev 1976; Paczynski 1992). In fact, the latter effect might even be responsible for explaining the highly super-Eddington nature (factor $\sim 100$) of accreting NSs in some ULXs; see Section 3.5. However, as mentioned in Section 4.8, the current spin distribution of observed DNS systems is fully consistent without the need for long-term accretion at rates exceeding ${\dot{M}}_{\mathrm{Edd}}$ by a large factor.

Applying a SN mass cut at the edge of the metal core and using the release of gravitational binding energy following Lattimer & Yahil (1989) give the mass of the newborn NS, ${M}_{\mathrm{NS}}={M}_{\mathrm{core},{\rm{f}}}-{M}_{\mathrm{NS}}^{\mathrm{bind}}$, where the released binding energy of the NS corresponds to a mass decrease of ${M}_{\mathrm{NS}}^{\mathrm{bind}}=0.084\,{M}_{\odot }\cdot {({M}_{\mathrm{NS}}/{M}_{\odot })}^{2}$. The resulting masses of the newborn NSs are in this specific case between 1.31 and $1.41\,{M}_{\odot }$ (see ${M}_{\mathrm{NS},2}$ in Table 5), given the relatively small initial helium star masses, ${M}_{\mathrm{He},{\rm{i}}}=3.0\,{M}_{\odot }$, assumed prior to Case BB RLO. Finally, knowing the mass of the newborn NS and all relevant pre-SN quantities, we are able to calculate the post-SN values of ${P}_{\mathrm{orb}}$.

At first look from Figure 8, one could be tempted to conclude that theoretical modeling with ${X}_{\mathrm{Edd}}\simeq 2.0$ can best reproduce the observational data. However, as we shall now demonstrate, the location of each observed DNS system in any $({P}_{\mathrm{orb}},P,e)$ plane is strongly dependent on both the mass of the progenitor of the exploding star (i.e., the helium star) and the magnitude and direction of the kick during the second SN.

To illustrate the spread in theoretically estimated points in the $({P}_{\mathrm{orb}},P)$ plane using different initial masses of the progenitor of the exploding star, we show in Figure 9 a number of points based on sample models from Tauris et al. (2015) for different initial helium star masses in the interval ${M}_{\mathrm{He},{\rm{i}}}=2.6\mbox{--}6.0\,{M}_{\odot }$ and with different initial orbital periods (see Table 6). In these models, the masses of the helium stars at the moment of the (electron capture or iron core collapse; see Section 6) SN explosions are $1.43\mbox{--}2.37\,{M}_{\odot }$, with metal cores of $1.37\mbox{--}2.15\,{M}_{\odot }$. The resulting NS masses cover a range between $1.2\mbox{--}1.9\,{M}_{\odot }$. This range is possibly wider by $\pm 0.1\,{M}_{\odot }$, if accounting for the uncertainty in the location of the mass cut and the NS EoS. For the theoretical values, shown as blue stars in Figure 9, we assumed a fixed value of ${X}_{\mathrm{Edd}}=2.0$. The fact that the scatter of the observed data points is smaller than that obtained from our applied range of theoretical models suggests some uniformity of the parent (pre-SN) systems. Note that following the method described above, we also find evidence for the formation of a few DNS systems outside the range of the plane shown in Figure 9, i.e., DNS systems with larger values of either ${P}_{\mathrm{orb}}$ or P.

Zoom In Zoom Out Reset image size

The presence of (even small) kicks also contributes significantly to smearing out the expected theoretical correlation between ${P}_{\mathrm{orb}}$ and P obtained for symmetric SNe. In the following sections, we illustrate the effect of an asymmetric SN.

In Figure 10, we have plotted the eccentricities of observed DNS systems (red points) as a function of their orbital periods. The measured error bars are too small to be seen. Also, a (weak) correlation between ${P}_{\mathrm{orb}}$ and eccentricity is expected from the binary stellar evolution during Case BB RLO (Tauris et al. 2015). Again, the argument relates to the decreasing amount of material transferred from helium star donors with increasing values of ${P}_{\mathrm{orb}}$ because of their shorter lifetimes prior to collapse—see Figure 16 in Tauris et al. (2015). As a result, these wide-orbit stars possess a large envelope mass by the time they explode which increases the instantaneous mass loss during the SN, thereby increasing the post-SN eccentricity (see the dashed blue line in Figure 10 as an illustrative example of symmetric SNe without a kick). However, this correlation also depends on the initial values of ${M}_{\mathrm{He},{\rm{i}}}$. Furthermore, it is more sensitive to even small kicks compared to the $({P}_{\mathrm{orb}},P)$ correlation.

Zoom In Zoom Out Reset image size

Indeed, the observed spread in eccentricities for a given value of ${P}_{\mathrm{orb}}$ can be understood if the second SN in these binaries (producing the second-born, and non-recycled, pulsar) is even slightly asymmetric, e.g., $w\lesssim 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$. In Figure 10, we simulated the dynamical effects of applying kicks of $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ to the stellar models presented in Table 5 with circular orbits prior to the SN. We performed 1000 trials with randomly oriented (isotropic) kicks. As can be seen from the figure, the resulting post-SN systems are spread out over a large area in the $({P}_{\mathrm{orb}},e)$ plane. An immediate result, however, is that with such small kicks, we cannot explain two of the DNS systems. These are the well-known PSR B1913+16 (the Hulse–Taylor pulsar; ${P}_{\mathrm{orb}}=0.323\,\mathrm{days}$ and e = 0.617) and PSR J1811−1736 (${P}_{\mathrm{orb}}=18.8\,\mathrm{days}$ and e = 0.828). As we shall discuss in more detail in Section 8, the latter DNS system can actually be reproduced from a small kick by increasing the value of ${M}_{\mathrm{He},{\rm{i}}}$, whereas PSR B1913+16 must have had $w\gtrsim 200\,\mathrm{km}\,{{\rm{s}}}^{-1}$ for any mass of the exploding star to match observations.

In relatively wide-orbit pre-SN binaries, even small kicks are able to produce a full range of post-SN eccentricities between $0\lt e\lt 1$, or even disrupt the system. This is demonstrated in Figure 11, where we modeled the explosions of ${M}_{\mathrm{He},{\rm{f}}}=3.0\,{M}_{\odot }$ stars with a fixed kick magnitude of $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and varied the direction of the applied kicks by systematically changing the two kick angles, θ and ϕ (see Section 6 for a definition). The pre-SN orbital period was assumed to be ${P}_{\mathrm{orb},{\rm{f}}}=25\,\mathrm{days}$ in all cases. It is seen that systems with θ less than some critical value, ${\theta }_{\mathrm{crit}}\simeq 83^\circ$, will always be disrupted—see discussions in Section 6. Another version of this plot, for an assumed ultra-stripped SN with a large kick in a tight orbit, is shown in the Appendix (Figure 40).

Zoom In Zoom Out Reset image size

In the previous sections, we have investigated the correlations for DNS systems between ${P}_{\mathrm{orb}}$ and P, and between ${P}_{\mathrm{orb}}$ and eccentricity, based on binary stellar evolution calculations starting from NS–helium star binaries. As a consequence, one might also expect a correlation between P and eccentricity. Indeed, such a correlation for DNS systems was suggested by McLaughlin et al. (2005) and Faulkner et al. (2005), and investigated theoretically by Dewi et al. (2005), based on the discovery of the first seven DNS systems, which could be fitted reasonably well onto a straight line in the $(P,e)$ plane.

Faulkner et al. (2005) suggested that in NS–helium star binaries, where a large amount of material is transferred from the helium star to the NS, leading to efficient spin up and thus a small value of P, less material will be ejected during the SN explosion, which results in smaller eccentricities. This qualitative argument is in good agreement with the more detailed analysis presented here and in Tauris et al. (2015).

In Figure 12, we show such an expected correlation between P and eccentricity for the same systems displayed in Figure 10 (see Table 5), and which subsequently exploded in either a symmetric SN (w = 0, dashed blue line) or in a SN with a fixed kick of $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The solid blue line represents the outcome of our theoretical models with applied kicks of $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ averaged over all kick directions that resulted in bound systems. We assumed again an Eddington accretion efficiency parameter ${X}_{\mathrm{Edd}}=2.0$ for the NS. It is clear that the observed systems (red points) do not follow the theoretical curve for a symmetric SN (w = 0). The reason is again, that even a small asymmetry in the SN will affect the eccentricities (see also the simulations and discussion in Dewi et al. 2005). This is also seen in Figure 12, where the black dots clustered on discrete vertical lines represent MC simulations for each of the 10 values of ${P}_{\mathrm{orb},{\rm{i}}}$, assuming a random (isotropic) kick direction for a fixed kick magnitude of $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (see also Figure 10). Note that considering only the eccentricities of the DNS systems as a function of P, we can explain (almost) all of the observed systems using $w\leqslant 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$. However, there are further observational constraints other than P, ${P}_{\mathrm{orb}}$, and eccentricity, such as proper motion, which in some cases require larger values of w (see discussions in Sections 6 and 8).

Zoom In Zoom Out Reset image size

In the previous four subsections (Sections 5.1–5.4), we have discussed the correlations between P, ${P}_{\mathrm{orb}}$, and eccentricity for DNS systems, as expected from the theoretical modeling of binary pre-SN stars and mainly under the assumption of purely symmetric SNe. However, we have also demonstrated (see Figures 9, 10, and 12) that applying different masses of the exploding star and/or taking into account (even small) kicks will camouflage any such correlation. We now continue discussing kicks in the context of a $(P,e)$ correlation.

5.4.1. Kick Magnitudes and Selection Effects

If the kick applied to the second-born NS is directed out of the orbital plane of the pre-SN system (i.e., if the kick angle $\phi \ne 0^\circ$ and $\phi \ne \pm 180^\circ$), then the spin axis of the recycled pulsar will be tilted with respect to the post-SN orbital angular momentum vector. This misalignment angle can be calculated as (Hills 1983)

$$\begin{eqnarray}&&\delta ={\tan }^{-1}\left(\displaystyle \frac{w\sin \theta \sin \phi }{\sqrt{{({v}_{\mathrm{rel}}+w\cos \theta )}^{2}+{(w\sin \theta \cos \phi )}^{2}}}\right)\end{eqnarray} \tag{ 16 }$$

and leads to the geodetic precession of the recycled pulsar (e.g., Kramer 1998). To use the observed misalignment angle to constrain the kick properties, we must rely on the assumption that accretion torques align the spin axis of the first-born NS with the orbital angular momentum vector during the recycling process (e.g., Hills 1983; Bhattacharya & van den Heuvel 1991). The fact that no observable effects of geodetic precession have been measured for pulsar A in the double pulsar system J0737−3039 implies that $\delta \lt 3\buildrel{\circ}\over{.} 2$ (Ferdman et al. 2013) and thus supports our assumption that ${\delta }_{i}=0$ prior to the second SN explosion in all DNS systems. Further observational evidence that such a spin–orbit alignment actually occurs during accretion in general was demonstrated for LMXBs by Guillemot & Tauris (2014), who found agreement between the viewing angles of binary MSPs (as inferred from γ-ray light-curve modeling) and their orbital inclination angles. Although the timescale of accretion during Case BB RLO in DNS progenitor systems is substantially shorter (by two to four orders of magnitude) than that in LMXBs, the torques at work will be larger due to the much higher mass-transfer rates during Case BB RLO (while the size of the magnetospheres remain roughly the same as a result of larger NS B-fields in DNS systems compared to fully recycled MSPs). Therefore, we find that it is reasonable to assume ${\delta }_{i}=0$ prior to the second SN explosion and thus legitimate to use the post-SN measurements of δ of the recycled NSs to constrain kicks in the second SN event.

Another important diagnostic quantity (aside from ${P}_{\mathrm{orb}}$, P, eccentricity, ${M}_{\mathrm{NS}}$, and δ) for understanding DNS formation is their post-SN systemic velocity, ${v}_{\mathrm{sys}}$. Any DNS system receives a recoil velocity relative to the center-of-mass rest frame of the pre-SN system. This is due to the combined effects of sudden mass loss and a kick velocity imparted on the newborn NS. From simple conservation of momentum considerations (e.g., following Tauris & Bailes 1996), we can write this 3D velocity as

$$\begin{eqnarray}&&{v}_{\mathrm{sys}}=\sqrt{{({\rm{\Delta }}{p}_{x})}^{2}+{({\rm{\Delta }}{p}_{y})}^{2}+{({\rm{\Delta }}{p}_{z})}^{2}}\,/\,({M}_{\mathrm{NS},1}+{M}_{\mathrm{NS},2}),\end{eqnarray} \tag{ 17 }$$

where the change in momentum is given by

$$\begin{eqnarray}\quad {\rm{\Delta }}{p}_{x} & = & {M}_{\mathrm{NS},2}\,w\cos \theta -{\rm{\Delta }}{{MM}}_{\mathrm{NS},1}\sqrt{G/({{Ma}}_{{\rm{i}}})},\\ \quad {\rm{\Delta }}{p}_{y} & = & {M}_{\mathrm{NS},2}\,w\sin \theta \cos \phi ,\\ \quad {\rm{\Delta }}{p}_{z} & = & {M}_{\mathrm{NS},2}\,w\sin \theta \sin \phi .\end{eqnarray} \tag{ 18 }$$

As an example, Figure 14 shows the calculated 3D systemic velocities from the simulations shown in Figure 10 where randomly oriented kicks with $w=50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ were applied. For these resulting systems, the average systemic velocity is $\langle {v}_{\mathrm{sys}}\rangle \simeq 30\,\mathrm{km}\,{{\rm{s}}}^{-1}$, which is still larger than the values observed for some of the DNS systems (Table 3) and thus indicating again that, at least in some cases, the kicks during the second SN have magnitudes of $w\lt 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (see detailed discussions further below).

Zoom In Zoom Out Reset image size

Finally, in addition to the kinematic effects from the second SN, DNS systems also have a relic systemic velocity component from the first SN. This value, however, is expected to be small. Based on the distribution and the peculiar velocities of HMXBs, Coleiro & Chaty (2013) find that their resulting systemic velocities are often only $\sim 10\mbox{--}20\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Even after taking into account the momentum absorbed by the massive companion star in HMXBs, we find that this result supports the hypothesis that exploding stars generally produce smaller kicks compared to the distribution of Hobbs et al. (2005), if they have lost (at least) their hydrogen-rich envelope via mass transfer prior to the SN (see Section 6.4 for discussions).

For purely symmetric SNe (w = 0), the equations governing the dynamical effects of the SN explosion (Equations (8)–(10)) simplify to (Flannery & van den Heuvel 1975)

$$\begin{eqnarray}&&\displaystyle \frac{{a}_{{\rm{f}}}}{{a}_{{\rm{i}}}}=\left[\displaystyle \frac{1-{\rm{\Delta }}M/M}{1-2{\rm{\Delta }}M/M}\right]\end{eqnarray} \tag{ 19 }$$

and

$$\begin{eqnarray}&&e={\rm{\Delta }}M/({M}_{\mathrm{NS},1}+{M}_{\mathrm{NS},2}),\end{eqnarray} \tag{ 20 }$$

and where the probability of remaining bound is always ${P}_{\mathrm{bound}}=1$ for ${\rm{\Delta }}M/M\lt 0.5$, whereas (following the virial theorem) all systems are disrupted if more than half of the total mass is lost, i.e., if ${\rm{\Delta }}M/M\gt 0.5$.

6.3.1. The Eccentricity Floor and the NS EoS

Whereas average kicks of $400\mbox{--}500\,\mathrm{km}\,{{\rm{s}}}^{-1}$ have been demonstrated for young isolated radio pulsars (Lyne & Lorimer 1994; Hobbs et al. 2005), it has been suggested for a couple of decades that exploding stars which are stripped in close binaries (i.e., SNe Ib/c) may produce substantially smaller NS kicks compared to SN II explosions of isolated, or very wide-orbit, stars (Tauris & Bailes 1996; Pfahl et al. 2002; Voss & Tauris 2003; Podsiadlowski et al. 2004; van den Heuvel 2004; Dewi et al. 2005; Bray & Eldridge 2016). This awareness came about from both theoretical and observational arguments, the former along the lines of stripped stars often leading to relatively fast explosions with small kicks (see Section 6.5), and the latter from comparison of the observed radio pulsar velocity distribution, or population synthesis studies based on this distribution, with the space velocities, orbital periods, and eccentricities of X-ray binaries, MSPs, and DNS systems. It is simply not possible to reproduce the observed data if exploding stars in close binaries in general would receive kicks of $400\mbox{--}500\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The conclusion also holds when selection effects are taken into account (Section 6.9).

Indeed, the evidence for small kicks has already been pointed out in several studies (e.g., Portegies Zwart & Yungelson 1998; Pfahl et al. 2002; Podsiadlowski et al. 2004; van den Heuvel 2004; Piran & Shaviv 2005; Schwab et al. 2010; Ferdman et al. 2013, 2014; Beniamini & Piran 2016). Moreover, for PSR J0737−3039 (the double pulsar) and PSR J1756−2251, Ferdman et al. (2013, 2014) derived small misalignment angles from observations of these systems, giving further support for small kicks.

On the other hand, we have clear evidence that even relatively large kicks can happen in close-orbit DNS progenitor systems. For example, to explain the characteristics (i.e., large proper motions) of the Hulse–Taylor pulsar (PSR B1913+16) and PSR B1534+12, a kick of at least $w\simeq 200\,\mathrm{km}\,{{\rm{s}}}^{-1}$ is needed (see Sections 8.4 and 8.12, and also, e.g., Wex et al. 2000; Wong et al. 2010). Furthermore, to explain the significant misalignment of the spin axis of the B-star companion from the orbital angular momentum vector in the pulsar binary system PSR J0045−7319, a large kick is needed, too (Kaspi et al. 1996; Wex 1998). Finally, there is some evidence for binaries being disrupted in the second SN, thereby explaining the observations of isolated mildly recycled radio pulsars with similar properties to the first-born NS in DNS systems (Lorimer et al. 2004). The velocities of such ejected NSs can be large even if the kick is small in case the former binary is tight and disrupted due to a large amount of mass loss during the SN event (Tauris & Takens 1998). However, ultra-stripping prior to the second SN event often prevents disruption due to mass loss, and in these cases a large kick is needed to break up the system.

Kicks are associated with explosion asymmetries and may arise from non-radial hydrodynamic instabilities in the collapsing stellar core (i.e., neutrino-driven convection and the standing accretion-shock instability; Blondin & Mezzacappa 2006; Scheck et al. 2006; Foglizzo et al. 2007; Marek & Janka 2009; Janka 2012, 2013; Bruenn et al. 2016; Foglizzo et al. 2015). These instabilities are thought to lead to large-scale anisotropies of the innermost SN ejecta, which interact gravitationally with the nascent NS and accelerate it on a timescale of several seconds (e.g., Janka 2012; Wongwathanarat et al. 2013).

The magnitude of the kick velocity imparted on a newborn NS is difficult to calculate, even for a given progenitor star modeled until the onset of core collapse. To predict the outcome of a stellar core collapse, O'Connor & Ott (2011) defined the compactness parameter of a stellar core at bounce as

$$\begin{eqnarray}&&{\xi }_{{\rm{M}}}=\displaystyle \frac{M/{M}_{\odot }}{R(M)/1000\,\mathrm{km}},\end{eqnarray} \tag{ 21 }$$

where R(M) is the radius of the coordinate enclosing mass, M. Sukhbold & Woosley (2014) argued that ξ at bounce is correlated with ξ at the onset of core collapse, and hence the calculated pre-SN density structures can help in forecasting the mapping between pre-SN (ideally even ZAMS for single stars) parameters and the outcome of a core collapse. Therefore, to enable an estimate of the expected kicks, it is necessary to calculate stellar density structures at the end of Case BB RLO evolutionary models (Tauris et al. 2013, 2015). So far, however, these progenitor models have not been calculated self-consistently to the very onset of core collapse (although see Moriya et al. 2017 for a composite calculation model integrating two stellar evolution codes).

In Wongwathanarat et al. (2013), a formula is presented (see their Equation (10)) for estimating maximum kick magnitudes. This formula is based on the assumption of a dipolar mass asymmetry, ${\rm{\Delta }}m$, between the two hemispheres, which is a fraction (maybe at most 20%–30%) of the total mass enclosed between the SN shock radius and the proto-NS at the time the explosion sets in. This mass can be determined numerically and scales with the mass in the postshock layer and thus, roughly, with the mass-infall rate of the collapsing stellar core as a progenitor-core dependent quantity. Moreover, the resulting kick depends on the shock velocity and thus on the explosion energy. It is important to stress that the actual magnitude of a considered NS kick is somewhere between zero and the maximum value derived from the formulae of Wongwathanarat et al. (2013), depending on the explosion asymmetry, which develops stochastically. The statistical probability distribution of this stochastic asymmetry must be determined using numerous explosion simulations. Finally, it has been suggested that multidimensional effects in the convective burning shells at the onset of core collapse can have an effect on the NS kick amplitude (e.g., Burrows & Hayes 1996; Arnett & Meakin 2011; Müller et al. 2016), i.e., pre-collapse asymmetries in the progenitor star could influence the asymmetry parameter of Equation (22) (see below) within the framework of the hydrodynamic kick mechanism.

Recently, Janka (2017) summarized the current theoretical understanding of how NS kicks come about using the gravitational tug-boat mechanism in asymmetric neutrino-driven core-collapse SN explosions. He derived a simple proportionality between the kick velocity, w, and the energy of the explosion, ${E}_{\exp }$:

$$\begin{eqnarray}&&w=211\,\mathrm{km}\,{{\rm{s}}}^{-1}\,C\left(\displaystyle \frac{{\alpha }_{\mathrm{ej}}}{0.1}\right)\left(\displaystyle \frac{{E}_{\exp }}{{10}^{51}\,\mathrm{erg}}\right){\left(\displaystyle \frac{{M}_{\mathrm{NS}}}{1.5{M}_{\odot }}\right)}^{-1},\end{eqnarray} \tag{ 22 }$$

where C is a constant typically of order unity, ${M}_{\mathrm{NS}}$ is the (baryonic) NS mass, and ${\alpha }_{\mathrm{ej}}$ is the momentum-asymmetry parameter of the explosion ejecta, whose statistics needs to be determined by hydrodynamic explosion modeling. The more asymmetric and more powerful the SN explosion is, the larger the NS kick can be (Figure 15). For example, according to Janka (2017) and references therein, core-collapse SNe of relatively massive iron cores have ${E}_{\exp }\simeq 3\times {10}^{50}\mbox{--}2.5\times {10}^{51}\,\mathrm{erg}$ and ${\alpha }_{\mathrm{ej}}=0\mbox{--}0.33$, yielding kick values in a broad range between $w=0\mbox{--}1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$, whereas progenitors of electron capture SNe (EC SNe) and low-mass iron core-collapse SNe (Fe CCSNe)—often relevant for ultra-stripped SNe, see below—typically have ${E}_{\exp }\simeq {10}^{50}\,\mathrm{erg}$ and ${\alpha }_{\mathrm{ej}}\lesssim 0.03$, resulting in very small kicks of $w\leqslant 10\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

We now continue discussing the kicks expected for ultra-stripped SNe. As we shall see in Section 8, the possibility of such small kicks is indeed in fine agreement with constraints from many observed DNS systems.

6.5.1. Ultra-stripped SN Kicks

While the binary stellar evolution arguments of Tauris et al. (2015) for small kicks hold best for the lowest-mass NSs formed in their models, it is possible, in principle, that even rather massive NSs are produced from ultra-stripped progenitors, provided that the whole (most of the) metal core collapses. Indeed, Tauris et al. (2015) predict that ultra-stripped SNe may potentially produce young NSs in DNS systems with a mass within the entire range $1.10\mbox{--}1.80\,{M}_{\odot }$. On the other hand, it is also expected that more massive pre-SN metal cores produce larger iron cores and thus more "normal" Fe CCSNe with larger explosion energies and therefore larger kicks.

6.6.1. Theoretical Arguments for a Mass-kick Relation

6.6.2. Observational Arguments for a Mass-kick Relation

In this regard, it is interesting to note that the DNS systems where the second-born NS has a relatively large mass $\gtrsim 1.33\,{M}_{\odot }$ are also those systems where it is evident that a large kick must have been at work; see Figure 16 and Table 7. The evidence for these large kick values is deduced from present post-SN orbital parameters of DNS systems; see Section 8. In Figure 16, we have excluded PSR J1518+4904 due to the relatively large error bar on its NS masses, and PSR J1906+0746 because of the lack of a constraint on the kick velocity (we find solutions for $w=0\mbox{--}1400\,\mathrm{km}\,{{\rm{s}}}^{-1}$, although with a peak at a low velocity of $w\simeq 70\,\mathrm{km}\,{{\rm{s}}}^{-1}$). For the second-born NSs with smaller masses, we note that the inferred kicks are quite small (or could potentially have been small).

Zoom In Zoom Out Reset image size

Table 7.  Young NS Masses and Resulting NS Kicks in DNS Systems

	${M}_{\mathrm{NS},2}$	${w}_{\mathrm{peak}}$	${w}_{\mathrm{interval}}$	e
Pulsar	(M⊙)	$(\mathrm{km}\,{{\rm{s}}}^{-1})$	$(\mathrm{km}\,{{\rm{s}}}^{-1})$
B1913+16	1.389	300	185−465	0.617
B1534+12	1.346	240	175−300	0.274
J1906+0746	1.291	70	0−1500	0.085
J0737−3039	1.249	18	0−70	0.088
J1756−2251	1.230	10	0−130	0.181
J0453+1559	1.174	45	15−175	0.113

Note. DNS systems listed in order of measured masses of the second-born (young) NSs. The values of ${w}_{\mathrm{peak}}$ (third column) denote the peak values of the solutions to their inferred kick velocity distribution with an interval of ${w}_{\mathrm{interval}}$ (fourth column), obtained from our analysis presented in Section 8. Note that the peak values were estimated using flat input distributions of pre-SN parameters and do not reflect weighted distributions based on binary stellar evolution. The orbital eccentricities are listed in the fifth column. See Table 2 for further characteristics.

Download table as:  ASCIITypeset image

Our hypothesis of a possible correlation between ${M}_{\mathrm{NS}}$ and w is still based on small number statistics. The number of observed DNS systems with precise NS mass measurements is limited. There could also be a transition region between small kicks and large kicks, possibly related to some element of stochasticity involved in the development of the explosion asymmetry (Wongwathanarat et al. 2013), which may mask any correlation.

To summarize our findings for NS kick magnitudes in the second SN in forming DNS systems, we conclude that while it is difficult to estimate kick magnitudes, all of the above arguments taken together allow us to speculate that most ultra-stripped SNe (EC SNe and, at least, Fe CCSNe with relatively small metal cores) generally lead to the formation of NSs with small kicks. Occasionally, however, large NS kicks are also at work in DNS systems.

To explain this, we advance the hypothesis that for Fe CCSNe, the kick magnitude may increase with the mass of the iron core of the exploding star, and thus the mass of the resulting NS. Since NS masses are affected by binary interactions of the progenitor system, kicks imparted on NSs in binaries are indeed expected to be different from those imparted on NSs born in isolation. Hence, if this hypothesis is correct, we can also use the inferred kick magnitudes from observations as a diagnostic tool for probing the evolutionary state of the stellar progenitor immediately prior to the second SN event.

Given our hypothesis of a mass-kick correlation for the second SN in DNS systems, one might also expect an (${M}_{\mathrm{NS}},e$) correlation, because systems that experience larger kicks are more likely to obtain larger eccentricities. However, as we have also demonstrated in this paper, different kick directions result in a large spread in post-SN eccentricities, even for relatively small kicks. Therefore, although the six DNS systems with precisely measured masses (Table 7) do show some hint of such a correlation (see plot in Figure 17), we would expect a rather large spread in the (${M}_{\mathrm{NS},2},e$) plane, as more DNS sources are discovered, albeit some average trend may persist. Observational selection effects related to eccentricities of DNS systems are discussed further in Section 9.2.

Zoom In Zoom Out Reset image size

There are two main reasons to believe that the majority of ultra-stripped SNe in DNS progenitor systems are Fe CCSNe rather than EC SNe. As discussed in detail in Podsiadlowski et al. (2004), Poelarends et al. (2008), Takahashi et al. (2013), Jones et al. (2013, 2014, 2016), Doherty et al. (2015), Woosley & Heger (2015), and Radice et al. (2017), the ZAMS mass range of isolated stars producing EC SNe is restricted to a narrow interval (possibly with a width of less than $1\,{M}_{\odot }$) located somewhere within the range $8\mbox{--}11\,{M}_{\odot }$. The corresponding mass window for producing EC SNe from helium stars in close binaries is restricted to a width of $\sim 0.2\,{M}_{\odot }$ (Tauris et al. 2015) and only for helium stars with a mass of $2.6\leqslant {M}_{\mathrm{He},{\rm{i}}}/{M}_{\odot }\leqslant 2.9$ (depending on the orbital period). Therefore, from a statistical point of view, since helium stars are formed in binaries with a broad range of masses, it would be somewhat unlikely that more than a few of the observed DNS systems experienced an EC SN (keeping in mind that both low-mass Fe CCSNe and EC SNe are expected to produce small kicks; see Section 6.5.1). Furthermore, it has been demonstrated (Tauris et al. 2013; Suwa et al. 2015; Moriya et al. 2017) that ultra-stripped helium donor stars with a final metal core mass just above the interval producing EC SNe ($\simeq 1.37\mbox{--}1.43\,{M}_{\odot }$) do actually produce Fe CCSNe. The initial helium star mass range for producing Fe CCSNe possibly extends up to $\sim 8\mbox{--}10\,{M}_{\odot }$ before BHs form, depending on, e.g., their wind mass-loss rates.

Newton et al. (2016) recently suggested that a relationship between the moment of inertia and the binding energy of NSs can be used in DNS systems formed via ultra-stripped SNe to distinguish between EC SNe and Fe CCSNe.

Finally, we notice that the observational classification of ultra-stripped SNe can be of both Type Ib and Type Ic (Tauris et al. 2015), depending on the amount of helium remaining in the envelope as well as the production and mixing of nickel during the explosion.

Binary systems that experience larger SN kicks may have a smaller probability of surviving the second SN. This provides a selection effect, favoring systems descending from small NS kicks to be overrepresented among the observed DNS systems. It is therefore natural to ask whether we can use the observed sample of DNS systems to probe the SN explosions.

In Figure 18, we have quantified this selection effect. It is seen that for tight pre-SN systems with ${P}_{\mathrm{orb}}=0.1\,\mathrm{days}$, the probability of surviving an ultra-stripped SN remains at 100% all the way up to $w\sim 220\,\mathrm{km}\,{{\rm{s}}}^{-1}$, whereas for a wide pre-SN system with ${P}_{\mathrm{orb}}=10\,\mathrm{days}$, the survival probability already declines from 100% at $w\sim 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$, and reaches zero, with no systems surviving, for $w\gt 330\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The survival probability for a more massive exploding star (${M}_{\mathrm{He},{\rm{f}}}=3.0\,{M}_{\odot }$, dotted line) is generally smaller because more mass is lost in the SN—this amount (${\rm{\Delta }}M$) is the sum of baryonic material ejected in the SN shell and the mass carried away by neutrinos due to the release of gravitational binding energy when the NS forms.

Zoom In Zoom Out Reset image size

Since the majority of the observed DNS systems seem to have had a pre-SN orbital period of $0.1\mbox{--}1\,\mathrm{days}$ (based on the observed distribution of ${P}_{\mathrm{orb}}$ for DNS systems and our simulations in Section 8), the survival probability remains within roughly 50% even for large kicks up to $400\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Therefore, we conclude that although there is indeed some selection at work for large kicks, the effect is rather limited. Moreover, we would expect a much broader distribution of eccentricities among the DNS population in case a large fraction of DNS systems almost become disrupted, in particular, for wide-orbit systems where gravitational damping has little effect.

The possibility that a fair fraction of bound DNS systems leave our Galaxy as a result of very large kicks also seems unlikely given that the local escape velocity in the Galactic rest frame is quite large, about $550\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Smith et al. 2007). Nevertheless, further population synthesis studies, e.g., of disrupted binaries (Lorimer et al. 2004), might help shed more light on the kick magnitudes.

For a measure of the post-SN merger timescale of a DNS system, ${\tau }_{\mathrm{gwr}}$ (Section 7), it is often useful to know the probability of the orbit shrinking as a consequence of the SN, i.e., $({a}_{{\rm{f}}}/{a}_{{\rm{i}}})\lt 1$, leading to small values of ${\tau }_{\mathrm{gwr}}$. In Figure 19, we have plotted these probabilities using Equation (13) for the same systems that exploded as in Figure 18. It is seen that the probability of obtaining a post-SN system with a reduced semimajor axis peaks at ∼34% (see below). For tight pre-SN systems, e.g., ${P}_{\mathrm{orb}}=0.1\,\mathrm{days}$, it is obvious that for small kick values ($w\lesssim 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$) the post-SN systems will always become wider (${P}_{{\rm{a}}}^{-}=0$). For very wide-orbit pre-SN systems, e.g., ${P}_{\mathrm{orb}}=10\,\mathrm{days}$, it is the large kicks (here $w\simeq 280\mbox{--}330\,\mathrm{km}\,{{\rm{s}}}^{-1}$) that always cause the post-SN systems to widen. In these cases, even if the kick is directed backwards toward the companion star, the kick "overshoots" the newborn NS into a wide orbit. For kicks above $330\,\mathrm{km}\,{{\rm{s}}}^{-1}$, the pre-SN systems with ${P}_{\mathrm{orb}}=10\,\mathrm{days}$ are all disrupted as shown in Figure 18.

Zoom In Zoom Out Reset image size

Comparing the solid, dashed, and dotted curves in Figure 19, we see that the more ultra-stripped the exploding star is, the more likely it is to reduce its orbit as a consequence of the SN explosion. The location and the values of the peak probabilities of ${P}_{{\rm{a}}}^{-}$ can easily be derived from Equation (13). We find that the location of the peak is at $(w/{v}_{\mathrm{rel}})=\sqrt{{\rm{\Delta }}M/M}$, and its value is given by

$$\begin{eqnarray}&&{P}_{{\rm{a}}}^{-}(\max )=\displaystyle \frac{1}{2}\left[1-\sqrt{\displaystyle \frac{{\rm{\Delta }}M}{M}}\right].\end{eqnarray} \tag{ 23 }$$

As an example, for the solid lines in Figure 19, we have $({M}_{\mathrm{He},{\rm{f}}}/{M}_{\odot },{\rm{\Delta }}M/{M}_{\odot })=(1.50,0.30)$ and ${M}_{\mathrm{NS},1}=1.40\,{M}_{\odot }$, which yields ${P}_{{\rm{a}}}^{-}(\max )=0.34$.

As we argued in Section 5.3, the wider the orbit of the exploding star the less stripped it will be prior to its collapse, and thus the larger ${\rm{\Delta }}M$. This leads to smaller probabilities for decreasing the orbit and therefore production of a DNS system that merges within a Hubble time. However, the situation is more complicated since the kick magnitude is likely to be correlated with ${\rm{\Delta }}M$, as argued in the previous subsections, which affects the probability of the post-SN orbit decreasing (Figure 19). Furthermore, the post-SN eccentricity is also important for ${\tau }_{\mathrm{gwr}}$ (see Section 7).

6.9.1. Selection Effects from the First SN Explosion

Finally, to complete the discussion of NS kicks and selection effects, and their implications for DNS production, we show in Figure 20 examples of the survival probabilities from the first SN explosion, i.e., the probability that the post-SN systems will remain bound (and later evolve to HMXB sources). The solid, dashed, and dotted lines are for different SNe Ib with typical values of $({M}_{\mathrm{He},{\rm{f}}}/{M}_{\odot },{M}_{2}/{M}_{\odot })$ expected for the first SN. It is seen that the probability of surviving the first SN is strongly dependent on the pre-SN orbital period (as is also the case for the second SN). This means that, from a pure kinematics point of view (disregarding the previous binary evolution), close-orbit HMXBs, like the RLO-HMXBs, are much more likely to form compared to wide-orbit HMXBs. In particular, we notice that there is a significant bias against the production of those HMXBs in orbits wide enough to survive the subsequent CE phase, which is needed to produce a DNS system (Section 3.4). For example, no binaries with pre-SN ${P}_{\mathrm{orb}}=120\,\mathrm{days}$ will survive a kick of $w\gt 300\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This selection effect must be kept in mind when trying to reconstruct the kick magnitudes based on observations of HMXB systems.

Zoom In Zoom Out Reset image size

For HMXBs to avoid merging in the CE phase, the systems must be in a fairly wide orbit prior to the HMXB phase. Based on a population synthesis study, Voss & Tauris (2003) concluded that this argument means that DNS systems can indeed only survive small kicks in the first SN in order to avoid disruption (indeed, Pfahl et al. 2002 discovered that there is a population of wide-orbit BeHMXBs with nearly circular orbits, indicating that in these systems the NS received small kicks). On the other hand, after the CE and spiral-in phases, the orbits of the surviving systems will always be tight, and therefore these systems will often be able to withstand even very large kicks in the second SN explosion.

A detailed population synthesis study to quantify the above-mentioned effects in more detail, and using weighted input distributions from binary stellar evolution, is currently ongoing (M. U. Kruckow et al. 2017, in preparation).

PSR J0453+1559 was recently discovered by Martinez et al. (2015). It is outstanding compared to the rest of the DNS population in that it has a large mass ratio between the two NS component masses of q = 1.33. The proper motion has been measured for this system, yielding a transverse velocity and thus a handle on its 3D systemic LSR velocity (Table 3). From our SN simulations, we are then able to place relatively tight constraints on, e.g., the kick magnitude of the second SN.

The simulated sample of solutions for this system is plotted in Figure 24, based on the NE2001 DM distance of 1.07 kpc. It is seen that there is a large probability for a relatively small imparted kick velocity of $w\lesssim 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ during the second SN (for an assumed upper limit on the resulting systemic velocity, ${v}_{\mathrm{sys}}\,\lt 85\,\mathrm{km}\,{{\rm{s}}}^{-1}$; Table 3), although the tail of solutions extends to relatively large kick values up to $w\simeq 175\,\mathrm{km}\,{{\rm{s}}}^{-1}$. For the mass of the exploding star (producing the second-born NS), we find solutions for ${M}_{\mathrm{He},{\rm{f}}}\lt 3.2\,{M}_{\odot }$. The pre-SN orbital period is found to be $2.5\lesssim {P}_{\mathrm{orb},{\rm{i}}}/\mathrm{days}\lesssim 4.5$.

Unfortunately, the distance to this pulsar remains rather uncertain, which makes it difficult to place tighter constraints on the kick magnitude. The kick is most likely to be $w\lt 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (although it could be somewhat larger). Using the Galactic electron distribution model of Yao et al. (2016) yields a DM distance of only 0.52 kpc. As a result, the systemic velocity of the system becomes smaller, and the imparted kick is found to be $w\lt 75\,\mathrm{km}\,{{\rm{s}}}^{-1}$. An upper limit for the mass of the exploding star is in this case ${M}_{\mathrm{He},{\rm{f}}}\lt 2.1\,{M}_{\odot }$. A parallax measurement would be able to resolve the distance discrepancy and thus allow for a tighter constraint on both w and ${M}_{\mathrm{He},{\rm{f}}}$. Any corrections for gravitational damping in this system are negligible due to its wide orbit (${P}_{\mathrm{orb}}=4.07\,\mathrm{days}$).

This DNS system is the famous double pulsar (Burgay et al. 2003; Lyne et al. 2004; Kramer et al. 2006b). It has been intensively timed since its discovery in 2003 and, due to its fortunate geometry (orbital inclination angle of 89°), its orbital and kinematic parameters are very well constrained. In particular, the misalignment angle of the first-born, recycled pulsar has been measured to be $\delta \lt 3\buildrel{\circ}\over{.} 2$ (Ferdman et al. 2013), and from VLBI observations (Deller et al. 2009) its distance ($\sim 1150\pm 200\,\mathrm{pc}$) and proper motion ($4.37\pm 0.65\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$) have been revealed. These measurements result in a small 2D transverse velocity of $\sim 25\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and, correcting for the unknown radial component and applying its location in the Galactic potential, we find an estimated 3D systemic LSR velocity of $\sim {33}_{-15}^{+24}\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (Table 3).

The resulting kick velocity found from our simulated solutions, where we restrict ourselves to only select systems with ${v}_{\mathrm{sys}}=18\mbox{--}57\,\mathrm{km}\,{{\rm{s}}}^{-1}$, is therefore $w\lesssim 70\,\mathrm{km}\,{{\rm{s}}}^{-1}$, as shown in Figure 25. As also seen in this figure, the kick magnitude distribution peaks at a value of only $\sim 18\,\mathrm{km}\,{{\rm{s}}}^{-1}$. A caveat to this conclusion is that the DNS system also received a (small) systemic velocity from the first SN, which was probably of the order of $10\mbox{--}20\,\mathrm{km}\,{{\rm{s}}}^{-1}$; see Section 6.2. A close inspection of Figure 25 shows that there are no solutions for ${v}_{\mathrm{sys}}\leqslant 32\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The reason for this is the finite amount of mass loss, ${\rm{\Delta }}M$.

Zoom In Zoom Out Reset image size

For a comparison, we plot in Figure 26 our simulated solutions without applying any constraints on ${v}_{\mathrm{sys}}$. The difference in the widths of the pre-SN solution parameter space is striking, especially for the value of the final mass of the exploding star, ${M}_{\mathrm{He},{\rm{f}}}$.

Zoom In Zoom Out Reset image size

Further important conclusions on the formation of PSR J0737−3039 can immediately be drawn from our simulations. Applying the constraint ${v}_{\mathrm{sys}}\lt 57\,\mathrm{km}\,{{\rm{s}}}^{-1}$ provides an almost unique solution to the pre-SN progenitor binary of PSR J0737−3039 (Figure 25), an argument that was first made by Piran & Shaviv (2005). This can be understood from the observations of the small proper motion of this system, which results in $(w/{v}_{\mathrm{rel}})\ll 1$, and thus a firm solution is obtained from the observed values of ${P}_{\mathrm{orb}}$ and e; see Equations (19) and (20). Similarly, it automatically follows that the misalignment angle in this case must be small (Equation (16)), in agreement with observations. For 94% of our solutions, the pre-SN binary had an orbital period of ${P}_{\mathrm{orb},{\rm{i}}}=0.085\pm 0.005\,\mathrm{days}$, and the mass of the (ultra-stripped) exploding star must have been ${M}_{\mathrm{He},{\rm{f}}}=1.56\pm 0.06\,{M}_{\odot }$. Interestingly enough, the very first calculation of a helium star–NS binary system leading to an ultra-stripped SN (Tauris et al. 2013) had pre-SN values of ${P}_{\mathrm{orb},{\rm{i}}}=0.070\,\mathrm{days}$ and ${M}_{\mathrm{He},{\rm{f}}}=1.50\,{M}_{\odot }$, and is thus basically a solution to the immediate progenitor of PSR J0737−3039. For this particular system, Tauris et al. (2013) found that the nature of the SN would be an Fe CCSN of Type Ic. From their binary stellar evolution modeling, we can therefore also deduce that the pre-SN system is the descendant of a helium star–NS binary with an initial (post-CE) orbital period of $\sim 0.1\,\mathrm{days}$ and a helium star ZAMS mass of $\sim 2.9\,{M}_{\odot }$.

In the discussion above, we have neglected any influence from gravitational damping on this system. The fact that we also observe the young pulsar in this system means that we can set an upper limit to the age of this binary. The characteristic spin-down age of pulsar B (the young NS) is about 50 Myr, which is also a typical radio pulsar lifetime for a non-recycled pulsar (Lorimer & Kramer 2004). Unless the braking index of this pulsar is significantly less than 3, this will serve as a proxy for the upper limit of its true age. However, in this particular binary, the magnetosphere of pulsar B, and thus its spin-down activity, is likely to be affected by the relativistic wind from pulsar A (the recycled pulsar, releasing spin-down energy at a rate of 3600 times more than pulsar B), as also indicated by observations (Lyutikov 2004; McLaughlin et al. 2004; Breton et al. 2012). Hence, the true age of pulsar B could possibly be more than 100 Myr (Lorimer et al. 2007).

For a DNS age of 100 Myr, we calculate corresponding post-SN parameters of ${P}_{\mathrm{orb}}=0.138\,\mathrm{days}$ and e = 0.119 (see location of red star), and the solution would be that the pre-SN binary had an orbital period of ${P}_{\mathrm{orb},{\rm{i}}}=0.11\pm 0.01\,\mathrm{days}$ and that the mass of the exploding star must have been ${M}_{\mathrm{He},{\rm{f}}}=1.59\pm 0.09\,{M}_{\odot }$. The distribution of valid kick velocities remains roughly similar to the one shown in Figure 25.

Since the exact age of PSR J0737−3039 remains unknown and could in principle take any value between 0 and 100 Myr (or even more), we conclude that this DNS system formed from a pre-SN binary with ${P}_{\mathrm{orb},{\rm{i}}}=0.10\pm 0.02\,\mathrm{days}$ and ${M}_{\mathrm{He},{\rm{f}}}=1.59\,\pm 0.09\,{M}_{\odot }$, as well as a small imparted kick, $w\lesssim 70\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

A final note on this double pulsar system is that it is often being used to constrain the empirical Galactic DNS merger rate. The outcome is dependent on the beaming geometry of both pulsars A and B. As an example, Kim et al. (2015) derived a value of ${21}_{-14}^{+28}\,{\mathrm{Myr}}^{-1}$.

Since the discovery of this DNS system (Nice et al. 1996), its component masses have remained difficult to constrain (Janssen et al. 2008). The total mass of the system ($M=2.72\,{M}_{\odot }$) remains a firm result, and we use that for our SN simulations of this system. Assuming that the mass of the second-born NS ${M}_{\mathrm{NS},2}$ falls within the range of such measured masses in other DNS systems ($1.17\mbox{--}1.39\,{M}_{\odot }$, Table 1), and using MC simulations with a flat probability distribution for ${M}_{\mathrm{NS},2}$ in this range, the mass of the first-formed NS is then simply evaluated from ${M}_{\mathrm{NS},1}=M-{M}_{\mathrm{NS},2}$. The results of our SN simulations are shown in Figure 27. The solutions yield values of the pre-SN orbital period between 4 and 12 days. The solutions for the mass of the exploding star in the second SN is restricted to $\lt 2.8\,{M}_{\odot }$, with a preference toward smaller values. Similarly, we only find solutions for a relatively small kick velocity ($w\lt 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$) with a peak²¹ near $w\simeq 40\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

Ongoing intensive data analysis, including data obtained in recent observation campaigns, yields preliminary values of ${M}_{\mathrm{NS},1}=1.41\,{M}_{\odot }$ and ${M}_{\mathrm{NS},2}=1.31\,{M}_{\odot }$ for the (mildly) recycled radio pulsar and its young NS companion, respectively (G. Janssen 2017, private communication). The reported error bars are still rather large, roughly $\pm 0.08\,{M}_{\odot }$. Redoing our MC simulations with these mass constraints does not introduce any significant changes to our results presented in Figure 27. Any corrections for gravitational damping in this system, during the time interval between its formation and now, are negligible due to its wide orbit (${P}_{\mathrm{orb}}=8.63\,\mathrm{days}$).

PSR B1534+12 was discovered by Wolszczan (1991) and studied further in detail by e.g., Stairs et al. (2002) and Fonseca et al. (2014). The system parameters are well-constrained (Table 2), not only because of its precisely measured masses, but also via its determined misalignment angle, $\delta =27\pm 3^\circ$, and its proper motion. The results from our SN simulations are shown in Figure 28. The mass of the exploding star is constrained to be ${M}_{\mathrm{He},{\rm{f}}}\lt 3.7\,{M}_{\odot }$. More noteworthy, the kick in the second SN is shown to have been large: $w\simeq 175\mbox{--}300\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This makes this binary special in the sense that it represents one out of only two DNS systems (besides PSR B1913+16) for which there is evidence for a large kick in an ultra-stripped SN.

Zoom In Zoom Out Reset image size

The reason that a large kick is required in this system is the relatively large measured value of δ. Without any constraints on this angle, solutions are found for a wide range of kick values, $50\lesssim w\lesssim 425\,\mathrm{km}\,{{\rm{s}}}^{-1}$, illustrating again the importance of measuring δ. The small error bar on the measurement of δ explains why the direction of this kick is seen to be confined to a relatively small area on the sphere surrounding the exploding star (lower-right panel). Corrections for gravitational damping in this system (since its formation) is rather limited due to the relatively small characteristic age of the recycled pulsar ($\tau =248\,\mathrm{Myr}$) and somewhat small eccentricity (e = 0.274); see the red solid star in the upper-left panel.

There are no mass constraints from observations of this DNS system (Keith et al. 2009). Assuming that the total mass of the system is within the interval of the total masses of all other known DNS systems ($M=2.57\mbox{--}2.88\,{M}_{\odot }$), and similarly that the mass of the second-born NS, ${M}_{\mathrm{NS},2}$, is also within the range ($1.17\mbox{--}1.39\,{M}_{\odot }$) measured in other DNS systems, we simulated the possible solutions for the pre-SN parameters and the NS kick. The results are shown in Figure 29. As can be seen, we cannot derive interesting constraints for the origin of this DNS system. The kick is likely to have been small ($\lt 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$), but we also find solutions for kicks up to $w=400\,\mathrm{km}\,{{\rm{s}}}^{-1}$. The wide orbit of this system excludes any gravitational damping of significance.

Zoom In Zoom Out Reset image size

This interesting DNS candidate was discovered by Ng et al. (2015). The nature of the system is still unclear. The spin period of the radio pulsar is 315 ms, which makes it quite likely to be a young pulsar. The young nature of this pulsar was recently confirmed by a measurement of $\dot{P}=2.43\times {10}^{-15}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ (C. Ng et al. 2017, in preparation). Indeed, if the NS had been even mildly recycled, we would have expected a smaller spin period, given its orbital period of 9.7 days (see Section 5.1). The compact object composition of the companion star remains unknown. It could be either a DNS system or a WD+NS system where the NS formed after the WD (Tauris & Sennels 2000). From the mass function, it is known that ${M}_{2}\gt 0.40\,{M}_{\odot }$, and given the unknown orbital inclination angle, the probability of a WD+NS binary is somewhat larger than that of a DNS system; see discussions in C. Ng et al. (2017, in preparation). On the other hand, the eccentricity of this system (e = 0.089) is similar to many DNS systems. The wide orbit of the system excludes any gravitational damping of significance, so this cannot be the explanation for its small eccentricity. Assuming that the system is indeed a DNS system, we show in Figure 30 our results from the simulations of the second SN. As in the case of PSR J1753−2240 discussed above, we have no constraints on the total mass nor on the value of ${M}_{\mathrm{NS},2}$, and we follow the same procedure as for PSR J1753−2240. The kick is likely to have been small ($\lt 50\,\mathrm{km}\,{{\rm{s}}}^{-1}$), but we also find solutions for kicks up to almost $w=400\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

PSR J1756−2251 (Faulkner et al. 2005) is very well constrained from precise measurements of its NS masses (${M}_{\mathrm{NS},1}=1.341\,{M}_{\odot }$, ${M}_{\mathrm{NS},2}=1.230\,{M}_{\odot }$) and the misalignment angle ($\delta \lt 34^\circ$) of the recycled pulsar (Ferdman et al. 2014). The proper motion measurement in declination, however, has been difficult to measure for this pulsar, since it is located near the ecliptic plane. Ferdman et al. (2014) derived a best estimate of ${\mu }_{\delta }=5.5\pm 7.0\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, which we combine with its proper motion in right ascension, ${\mu }_{\alpha }=-2.42\pm 0.08\,\mathrm{mas}\,{\mathrm{yr}}^{-1}$, to estimate its transverse velocity and a determination of its 3D systemic velocity of ${v}^{\mathrm{LSR}}=41\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (${v}_{1\,\sigma }^{\mathrm{LSR}}=21\mbox{--}71\,\mathrm{km}\,{{\rm{s}}}^{-1}$; see Table 3).

This DNS system shares many properties with the double pulsar (PSR J0737−3039) and can be seen as a wider-orbit version of it. The SN simulations reveal a similar solution for the origin of the second-born NS as an ultra-stripped exploding star with a mass of $1.72\pm 0.22\,{M}_{\odot }$. The kick velocity is found to be $w\lt 130\,\mathrm{km}\,{{\rm{s}}}^{-1}$, with a main peak in solutions toward $w\to 0$. It is interesting to note that this DNS system is the only one for which half of the solutions require a kick in a forward direction, i.e., $\theta \lt 90^\circ$. This is partly due to the relatively small 3D systemic velocity inferred for of this system.

Using the NE2001 electron distribution model yields a DM distance of $d=0.73\,\mathrm{kpc}$, whereas applying the model of Yao et al. (2016) results in a distance four times larger ($d=2.81\,\mathrm{kpc}$). Such a large distance would significantly increase our derived interval of possible kick velocities for this system. However, a firm upper value for the distance of 1.2 kpc can be derived using the measured orbital period decay (Ferdman et al. 2014). This example illustrates again the need for precise distance determinations via parallax measurements.

The orbital period of this system of ${P}_{\mathrm{orb}}=0.320\,\mathrm{days}$, and its eccentricity of e = 0.181, means that some gravitational damping may have been at work since its formation. Assuming an age of $t=\tau =443\,\mathrm{Myr}$ yields initial values of ${P}_{\mathrm{orb}}=0.353\,\mathrm{days}$ and e = 0.198 right after the second SN (see the red solid star in the upper-left panel in Figure 31). However, this potential slight change in initial values does not impose any effects on the constraints of the pre-SN parameters and the kick velocity derived here.

PSR J1811−1736 (Corongiu et al. 2007) is a wide-orbit, highly eccentric DNS system with ${P}_{\mathrm{orb}}=18.8\,\mathrm{days}$ and e = 0.828. Its total mass is $2.57\,{M}_{\odot }$, whereas its individual NS mass components remain unknown. Assuming that the mass of the second-born NS, ${M}_{\mathrm{NS},2,}$ is within the range of observed masses in other DNS systems ($1.17\mbox{--}1.39\,{M}_{\odot }$), we simulated the possible solutions for the pre-SN parameters and the NS kick. The results are shown in Figure 32. As can be seen, there is a large range of possible values for all pre-SN parameters and the kick, and we conclude that we cannot derive interesting constraints for the origin of this DNS system. The wide orbit of PSR J1811−1736 excludes any gravitational damping of significance.

Zoom In Zoom Out Reset image size

This DNS system was discovered by Champion et al. (2004). The total mass of the system is $2.59\,{M}_{\odot }$ and the second-born NS has a mass ${M}_{\mathrm{NS},2}\gt 1.22\,{M}_{\odot }$. Assuming as usual ${M}_{\mathrm{NS},2}\leqslant 1.39\,{M}_{\odot }$, our SN simulations yield the results shown in Figure 33. Also in this case, there is a large range of possible values for all pre-SN parameters and the kick, and thus we cannot derive interesting constraints for the origin of this DNS system. The relatively wide orbit (${P}_{\mathrm{orb}}=1.18\,\mathrm{days}$), and the small eccentricity (e = 0.139), of the system means that its orbital parameters have hardly changed since its formation.

Zoom In Zoom Out Reset image size

This is a remarkable binary since the observed radio pulsar is the young, second-born NS (Lorimer et al. 2006). The evidence is a combination of a large B-field (i.e., $\dot{P}=2.03\times {10}^{-14}\,{\rm{s}}\,{{\rm{s}}}^{-1}$) and a relatively slow spin period ($P=144\,\mathrm{ms}$); see Table 2 and Figure 2. These values are typical for a young, non-recycled NS. The NS component masses are measured to be ${M}_{\mathrm{NS},1}=1.322\,{M}_{\odot }$ and ${M}_{\mathrm{NS},2}=1.291\,{M}_{\odot }$, respectively (van Leeuwen et al. 2015). Since the companion star remains unseen, this system could, in principle, also be a WD+NS system. ONeMg WD masses of $1.322\,{M}_{\odot }$ are certainly possible. We show our simulations of the (second) SN in Figure 34. However, we cannot make any interesting constraints on the pre-SN parameters of the system. A remarkable result though is that this binary could, in principle, have survived a NS kick velocity of about $w=1500\,\mathrm{km}\,{{\rm{s}}}^{-1}$. A proper motion measurement would be able to constrain the kick velocity. Since we treat the first-formed compact object as a point mass in these simulations, the results are equally valid for a WD+NS system where the NS forms after the WD (Tauris & Sennels 2000).

Zoom In Zoom Out Reset image size

The close orbit of this system (${P}_{\mathrm{orb}}=0.166\,\mathrm{days}$) implies that gravitational damping since birth could be of significance. However, this NS must be quite young given its high value of $\dot{P}=2.03\times {10}^{-14}\,{\rm{s}}\,{{\rm{s}}}^{-1}$ ($\tau =114\,\mathrm{kyr}$), and therefore its present orbital parameters resembles very well its birth parameters. O'Shaughnessy & Kim (2010) argued that this DNS system makes quite a significant contribution to the empirical estimate of the Galactic DNS merger rate.

Lazarus et al. (2016) recently announced the discovery of this DNS system with ${P}_{\mathrm{orb}}=0.206\,\mathrm{days}$ and e = 0.090. The total mass of the system is reported to be $M=2.88\,{M}_{\odot }$. The NS masses are expected to be measured soon, but until then we assume that the mass of the second-born NS, ${M}_{\mathrm{NS},2},$ is within the usual range of 1.17–1.39 M_⊙, as measured in other DNS systems. In Figure 35, we show our simulated possible solutions for the pre-SN parameters and the NS kick. As can be seen, we cannot derive interesting constraints for the origin of this DNS system. The orbit of the binary could have been affected significantly by gravitational damping since its formation. The main reason for this is the large characteristic age of the system ($\tau =2690\,\mathrm{Myr}$), which potentially allows for a large orbital decay on a gigayear timescale. In the upper-left panel, we have plotted its birth location with a red solid star (${P}_{\mathrm{orb}}=0.438\,\mathrm{days}$ and e = 0.189) for an assumed age $t=\tau$. These initial values are almost independent of the exact NS masses (here we simply assumed ${M}_{\mathrm{NS},1}=1.58\,{M}_{\odot }$ and ${M}_{\mathrm{NS},2}=1.30\,{M}_{\odot }$, motivated by typical values of ${M}_{\mathrm{NS},2}$; see Table 1).

Zoom In Zoom Out Reset image size

The Hulse–Taylor pulsar (Hulse & Taylor 1975; Kramer 1998; Weisberg et al. 2010) is an unusual case since this DNS system must have been accompanied by a rather large kick during the second SN explosion (e.g., Cordes & Wasserman 1984; Burrows & Woosley 1986; Wex et al. 2000; Willems et al. 2004). We confirm this result from our analysis of its 3D systemic velocity (Section 2.2) and the SN simulations presented in Figure 36. We find kick solutions for the entire interval of $w=185\mbox{--}465\,\mathrm{km}\,{{\rm{s}}}^{-1}$. This result (as well as that for the other DNS system with a large kick, PSR B1534+12) is in excellent agreement with the findings of Wong et al. (2010). For the mass of the exploding star, however, we cannot constrain its mass to better than ${M}_{\mathrm{He},{\rm{f}}}\lt 4.4\,{M}_{\odot }$. Indeed, the possibility of a somewhat more massive exploding star, compared to some of the DNS systems that experienced small kicks (e.g., PSR J0737−3039 and PSR J1756−2251), may be a clue to understanding the large kick in this system (see discussions in Section 6.6).

Zoom In Zoom Out Reset image size

A caveat to this conclusion is that the distance to PSR B1913+16 is uncertain. The Yao et al. (2016) DM distance of 5.25 kpc is only half the distance of 9.8 kpc provided by Weisberg et al. (2008) and Weisberg & Huang (2016), based on the Galactic model of Reid et al. (2014). Adapting this smaller distance results in typical kick velocities that are about $50\mbox{--}100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ smaller than those quoted above, although with a minimum kick velocity remaining at a relatively large value of $w=165\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Perhaps more interesting and important, the mass of the exploding star can in this case be constrained to be ${M}_{\mathrm{He},{\rm{f}}}\lt 2.6\,{M}_{\odot }$.

The torus of possible kick angles in the (ϕ, θ) plane is a result of the constraint from the measured misalignment angle, $\delta =18\pm 6^\circ$ (Ferdman et al. 2013). The combination of a small orbital period of ${P}_{\mathrm{orb}}=0.323\,\mathrm{days}$ and a large eccentricity of e = 0.617 means that gravitational damping is significant in this system, as discussed in Section 7. We illustrate this by marking the birth location of PSR B1913+16 in the (${P}_{\mathrm{orb}},e$) plane (upper-left panel in Figure 36) with a red solid star for a DNS age of $t=108.7\,\mathrm{Myr}$. The SN simulations leading to these birth properties are shown in Figure 37. The differences in the solutions for the pre-SN parameters and the second NS kick are very limited compared to the solutions shown in Figure 36. For an age of $t=108.7\,\mathrm{Myr}$, we do find solutions for a slightly more massive exploding star (up to $4.8\,{M}_{\odot }$). However, the most likely kick remains at roughly $200\mbox{--}450\,\mathrm{km}\,{{\rm{s}}}^{-1}$.

Zoom In Zoom Out Reset image size

This important very wide-orbit DNS system was discovered by Swiggum et al. (2015). The binary has ${P}_{\mathrm{orb}}=45\,\mathrm{days}$ and a marginally recycled pulsar spinning at $P=185\,\mathrm{ms}$. The total mass of this system is $2.59\,{M}_{\odot }$ and the second-born NS has a mass ${M}_{\mathrm{NS},2}\gt 1.30\,{M}_{\odot }$. Assuming as usual ${M}_{\mathrm{NS},2}\leqslant 1.39\,{M}_{\odot }$, our SN simulations yield the results shown in Figure 38. Also in this case, there is a large range of possible values for all pre-SN parameters and the kick, and thus we cannot derive interesting constraints for the origin of the DNS system. The kick was most likely rather small ($w\lt 100\,\mathrm{km}\,{{\rm{s}}}^{-1}$), although we find a small tail of possible solutions up to $290\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Gravitational damping is not relevant for this system. Assuming an age of 50 Myr or 100 Myr reveals an initial spin period of the marginally recycled pulsar of 154 ms or 115 ms, respectively, assuming a spin evolution with a braking index of n = 3 (see Equation (29)). The orbital parameters and the spin of this marginally recycled pulsar are discussed further in Tauris et al. (2015).

Zoom In Zoom Out Reset image size

The extreme stripping of helium stars during Case BB RLO, following the CE phase, leads to relatively low masses of the second exploding star in DNS systems—typically of values ${M}_{\mathrm{He},{\rm{f}}}\leqslant 2.5\,{M}_{\odot }$ for a broad range of initial values of orbital periods and helium star masses (Tauris et al. 2015). The solutions for the second SNe presented in Figures 24–38 are based purely on kinematics. Many of the solutions with massive exploding stars (up to $\sim 7\,{M}_{\odot }$) are therefore not physical, since binary stellar evolution would, in most cases, prevent such massive exploding stars. Only in a few cases is it possible that very massive helium stars (WR stars) may produce exploding stars with such large masses if their stellar wind mass-loss rate is small (e.g., Yoon et al. 2010 and references therein). To illustrate the strong dependence on stellar wind mass-loss rate, Tauris et al. (2015) calculated that the outcome of a binary system with a NS and a $10\,{M}_{\odot }$ helium star could be an exploding star with any mass between $3\mbox{--}7\,{M}_{\odot }$. However, if such a binary system is in a close orbit, GW radiation will force the system into RLO while the helium star is still massive, causing the system to merge upon RLO. Alternatively, for a wide-orbit system, the binary may expand significantly as a result of the strong wind such that the NS escapes being recycled via Case BB RLO.

Since some of the quantities plotted in Figures 24–38 are correlated with the mass of the exploding star (e.g., the magnitude and required direction of the kick velocities leading to a given DNS solution), care must be taken when drawing conclusions from these SN simulations. Full population synthesis is needed to disentangle cross-dependent quantities and enable more trustworthy probability distributions of various parameters at play.

In Figure 23, we have plotted the directions of all kick velocities leading to a solution for the simulations of the second SN of the DNS systems shown in Figures 24–38. We restricted the number of solutions to 3000 for each system and stacked together the ($\phi ,\theta$) solution planes on top of each other. The darker the blue-scale color (log scale), the more solutions in a given direction. A remarkable apparent preference is seen for a kick direction out of the plane of the pre-SN binary. The pre-SN orbital plane is along $\phi =0^\circ$ and $\phi =\pm 180^\circ$, and the kick solutions clearly cluster in directions out of the plane ($\phi \simeq 90^\circ$ and $\phi \simeq -90^\circ$). Furthermore, the kicks leading to the solutions show a significant preference for a direction backwards ($\theta \gt 90^\circ$). This latter result is perhaps not surprising since the probability of the system remaining bound is larger for a backward kick; see Section 6. For a verification of an isotropic input (trial) distribution of kick velocity directions, see Figure 39 in the Appendix.

Zoom In Zoom Out Reset image size

Given the above discussion on kinematic solutions versus realistic (from a binary stellar evolution point of view, Section 9.1) masses of the exploding stars (${M}_{\mathrm{He},{\rm{f}}}$), we performed an additional simulation of kick distributions in the observed DNS systems, in which we only kept solutions with ${M}_{\mathrm{He},{\rm{f}}}\lt 2.5\,{M}_{\odot }$. The outcome of this experiment remains fairly identical to the result shown in Figure 23. The only minor difference is that the peak of the kick directions is located at a slightly smaller value of $\theta \simeq 125^\circ$, and that more solutions allow for kicks in a forward direction ($\theta \lt 90^\circ$). The latter result is expected for smaller masses of the exploding stars.

Further motivated by the puzzle of an apparent NS kick anisotropy, we performed additional investigations and found that this anisotropy can be understood from a combination of an input prior and a selection effect.

9.2.1. Kick Anisotropies and Input Priors

First, we discuss how our application of unconstrained input distributions of kick magnitudes (i.e., without upper limits) can affect the solutions of kick angles for those DNS systems in which there are no observational constraints on the proper motion. These DNS systems include PSRs J1753−2240, J1755−2550, J1829+2456, J1906+0746, J1913+1102, and J1930−1852, which have small ($e\simeq 0.1$) or intermediate ($e\simeq 0.3\mbox{--}0.4$) eccentricities, and for which our SN simulations revealed kick solutions up to several $100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (in a few cases, even above $1000\,\mathrm{km}\,{{\rm{s}}}^{-1}$). In order for the simulated progenitor systems to remain bound and reproduce the observed DNS systems, the kick directions must preferentially be in the regions near the yellow crosses in Figure 23. The reason for this is that to compensate for the large applied kick magnitudes, the kick directions are limited in order to reproduce the small post-SN eccentricities (see Figures 11 and 40), which are needed to match most of the observed DNS systems. If instead we impose limitations on the maximum kick magnitude in the second SN in these systems, i.e., treating the kick magnitudes as elicited priors, and only allow kicks up to some limited maximum, say $w=50\mbox{--}60\,\mathrm{km}\,{{\rm{s}}}^{-1}$, then the apparent kick anisotropy along ϕ disappears.

Zoom In Zoom Out Reset image size

If assuming that nature produces an isotropic distribution of kick directions, we can therefore reverse the above argument and claim further evidence for small kicks in ultra-stripped SNe in general (but not always, as demonstrated by PSRs B1913+16 and B1534+12).

9.2.2. Kick Anisotropies and Selection Effects

Merging DNS and NS+BH binaries are considered the most promising sources to explain the phenomenology of short GRBs (Paczynski 1986; Eichler et al. 1989; Narayan et al. 1992; Fox et al. 2005; Rezzolla et al. 2011; Berger 2014). These short GRBs (sGRBs) are possibly followed by a so-called "kilonova" or "macronova" (Li & Paczyński 1998; Kulkarni 2005; Metzger et al. 2010; Barnes & Kasen 2013; Berger et al. 2013; Tanaka & Hotokezaka 2013; Tanvir et al. 2013; Fernández & Metzger 2016; Metzger 2017; Gottlieb et al. 2017). According to the standard scenario for the generation of sGRBs, an accretion-powered jet from a remnant BH torus system is formed within $10\mbox{--}100\,\ \mathrm{ms}$ after the merger event (e.g., Narayan et al. 1992; Janka et al. 1999; Rezzolla et al. 2011; Just et al. 2016). Moreover, the energy release should cease once the torus has been accreted on a timescale of $\leqslant 1\,{\rm{s}}$, which is consistent with the duration of the sGRB prompt emission of less than $\sim 2\,{\rm{s}}$.

Observations by the Swift satellite, however, have recently revealed long-lasting (${10}^{2}\mbox{--}{10}^{5}\,{\rm{s}}$) X-ray activity in a large fraction of sGRB events, which is indicative of a continuous energy release from a central engine on timescales orders of magnitude larger than the typical torus accretion timescale (Rowlinson et al. 2010, 2013; Gompertz et al. 2013, 2014; Lü et al. 2015). This late X-ray activity is difficult to explain with the accretion of the torus or the prolonged afterglow radiation produced by the interaction of the jet with the ambient medium (Kumar & Zhang 2015). Rezzolla & Kumar (2015) and Siegel & Ciolfi (2016) recently suggested a solution to this issue. If a large fraction of DNS mergers result in the formation of a metastable or long-lived NS rather than a BH, extraction of rotational energy via magnetic spin-down from such an object (typically a millisecond magnetar; Usov 1992) could power the observed long-lasting X-ray activity. Such a model will not work for the NS+BH scenario for sGRBs and would therefore exclude a large fraction of sGRBs from originating from NS+BH mergers.

The merger process of DNS systems is a vital source of chemical enrichment of the interstellar medium with heavy r-process elements such as, for example, gold (e.g., Eichler et al. 1989; Freiburghaus et al. 1999; Wanajo et al. 2014; Just et al. 2015; Rosswog 2015; Hirai et al. 2015; Beniamini et al. 2016b; Wu et al. 2016). These elements were previously thought to be synthesized primarily in SNe. The total mass of the ejected baryonic matter from DNS mergers, including the yield of heavy r-process elements, is typically of the order $0.01\,{M}_{\odot }$, depending on the mass ratio between the two merging NSs. Until recently, all observed DNS systems had a mass ratio $q\simeq 1$. The discovery of PSR J0453+1559 (Martinez et al. 2015) revealed the first known DNS system with a more asymmetric pair of NS masses (1.56 and $1.17\,{M}_{\odot }$) and $q\gt 1.3$. This mass asymmetry is very important as it may possibly increase (see below) the yield of expelled heavy r-process elements, when the lighter NS fills its Roche lobe and becomes tidally disrupted.

In Tauris et al. (2015), the modeling of DNS progenitor systems resulted in second-born NS masses in the interval $1.1\mbox{--}1.8\,{M}_{\odot }\,(\pm 0.1\,{M}_{\odot })$. The fact that the lightest observed first-born NSs in DNS systems have masses $\lt 1.2\,{M}_{\odot }$ means that, in principle, we might even expect a future detection of a DNS merger with a mass ratio $q\geqslant 1.5$. Such a DNS merger would produce r-process rich ejecta with larger asymmetry and a different direction-dependent composition than those expected for a DNS merger with $q\simeq 1$ (Just et al. 2015), and possibly, but not necessarily, a larger yield (depending on the NS masses and the NS EoS; Just et al. 2015; Lehner et al. 2016).

The Galactic merger rate of DNS systems has been estimated for almost four decades (e.g., Clark et al. 1979). As we have discussed in the previous sections, the standard formation scenario of DNS binaries involves a number of highly uncertain aspects of binary interactions (e.g., Voss & Tauris 2003; Tauris & van den Heuvel 2006; Belczynski et al. 2008; Dominik et al. 2012). The main uncertainties include, in particular, the treatment of CE evolution (Ivanova et al. 2013; Kruckow et al. 2016), SN kicks (Janka 2012, 2017) and the efficiency of accretion and spin-up from mass transfer (e.g., Tauris et al. 2015 and this paper). Together, these processes lead to large uncertainties in the expected merger rates of several orders of magnitude. Currently, the DNS merger rate is predicted to be about $1\mbox{--}100\,{\mathrm{Myr}}^{-1}$ per Milky Way Equivalent Galaxy (Abadie et al. 2010).

The vast majority,²² almost all, of DNS systems that merge within a Hubble time due to GW radiation are produced by ultra-stripped SNe (Tauris et al. 2015). And because some relatively wide-orbit DNS systems are also produced via ultra-stripped SNe, the observed rate of ultra-stripped SNe is expected to be larger than the rate of DNS mergers. Furthermore, ultra-stripped SNe may also originate in binary systems with a BH or a WD accretor, and thus the total rate of ultra-stripped SNe is estimated to be one in every 100–1000 SNe (Tauris et al. 2013).

After the recent success in detecting GW signals from merging BHs by advanced LIGO (Abbott et al. 2016a, 2016b, 2017), it is also expected that GWs from merging DNS systems will be detected in the near future with LIGO and other detectors like VIRGO and KAGRA. Thus, it will soon be possible to obtain DNS merger rate constraints from GW observations. This should place a lower limit on the rate of ultra-stripped SNe. Comparing the rates of ultra-stripped SNe inferred from high-cadence optical surveys and DNS mergers from GW detectors, we will be able to test whether the ultra-stripped SN channel is actually the major path for forming DNS systems, and thus constrain binary stellar evolution using GWs.