# Formation of Double Neutron Star Systems

T.M. Tauris, M. Kramer, P.C.C. Freire, N. Wex, H.-T Janka, N. Langer, Ph.

Podsiadlowski, E. Bozzo, S. Chaty, M.U. Kruckow, et al.

► To cite this version:

T.M. Tauris, M. Kramer, P.C.C. Freire, N. Wex, H.-T Janka, et al. Formation of Double Neutron Star Systems. The Astrophysical Journal, 2017, 846 (2), pp.170. 10.3847/1538-4357/aa7e89. hal-01645592

HAL Id: hal-01645592

https://hal.science/hal-01645592v1

Submitted on 10 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Formation of Double Neutron Star Systems

T. M. Tauris<sup>1,2</sup>①, M. Kramer<sup>1</sup>①, P. C. C. Freire<sup>1</sup>①, N. Wex<sup>1</sup>, H.-T. Janka<sup>3</sup>①, N. Langer<sup>2</sup>, Ph. Podsiadlowski<sup>2,4</sup>, E. Bozzo<sup>5</sup>①, S. Chaty<sup>6,7</sup>①, M. U. Kruckow<sup>2</sup>, E. P. J. van den Heuvel<sup>8</sup>, J. Antoniadis<sup>1,9</sup>①, R. P. Breton<sup>10</sup>①, and D. J. Champion<sup>1</sup>①

<sup>1</sup>Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, D-53121 Bonn, Germany; tauris@astro.uni-bonn.de

<sup>2</sup>Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany

<sup>3</sup>Max Planck Institut für Astrophysik, Karl-Schwarzschild-Strasse 1, D-85748 Garching, Germany

<sup>4</sup>Department of Astronomy, Oxford University, Oxford OX1 3RH, UK

<sup>5</sup>ISDC Data Centre for Astrophysics, Chemin d'Ecogia 16, CH-1290 Versoix, Switzerland

<sup>6</sup>Laboratoire AIM (UMR 7158 CEA/DSM-CNRS-Université Paris Diderot), Irfu/Service d'Astrophysique, Centre de Saclay, F-91191 Gif-sur-Yvette Cedex, France

<sup>7</sup>Institut Universitaire de France, 103 Boulevard Saint-Michel, F-75005 Paris, France

<sup>8</sup>Astronomical Institute Anton Pannekoek, University of Amsterdam, P.O. Box 94249, 1090 GE Amsterdam, The Netherlands

<sup>9</sup>Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada

<sup>10</sup>Jodrell Bank Centre for Astrophysics, School of Physics and Astronomy, The University of Manchester, Manchester M13 9PL, UK

Received 2017 March 20; revised 2017 June 9; accepted 2017 June 21; published 2017 September 13

## Abstract

Double neutron star (DNS) systems represent extreme physical objects and the endpoint of an exotic journey of stellar evolution and binary interactions. Large numbers of DNS systems and their mergers are anticipated to be discovered using the Square Kilometre Array searching for radio pulsars, and the high-frequency gravitational wave detectors (LIGO/VIRGO), respectively. Here we discuss all key properties of DNS systems, as well as selection effects, and combine the latest observational data with new theoretical progress on various physical processes with the aim of advancing our knowledge on their formation. We examine key interactions of their progenitor systems and evaluate their accretion history during the high-mass X-ray binary stage, the common envelope phase, and the subsequent Case BB mass transfer, and argue that the first-formed NSs have accreted at most ~0.02 $M_{\odot}$. We investigate DNS masses, spins, and velocities, and in particular correlations between spin period, orbital period, and eccentricity. Numerous Monte Carlo simulations of the second supernova (SN) events are performed to extrapolate pre-SN stellar properties and probe the explosions. All known close-orbit DNS systems are consistent with ultra-stripped exploding stars. Although their resulting NS kicks are often small, we demonstrate a large spread in kick magnitudes that may, in general, depend on the past interaction history of the exploding star and thus correlate with the NS mass. We analyze and discuss NS kick directions based on our SN simulations. Finally, we discuss the terminal evolution of close-orbit DNS systems until they merge and possibly produce a short γ-ray burst.

**Key words:** gravitational waves – pulsars: general – stars: mass-loss – stars: neutron – supernovae: general – X-rays: binaries

## 1. Introduction

The evolution of massive binary stars and the subsequent production of pairs of compact objects in tight orbits play a central role in many areas of modern astrophysics, including the origin of different types of supernova (SN) explosions (Whelan & Iben 1973; Iben & Tutukov 1984; Yoon et al. 2010), possibly including long gamma-ray bursts (GRBs; Woosley 1993; Cantiello et al. 2007), the modeling of accretion processes in X-ray binaries (Lewin & van der Klis 2006), and the formation of millisecond pulsars (MSPs; Bhattacharya & van den Heuvel 1991; Tauris & van den Heuvel 2006; Patruno & Watts 2012). In some cases, we even expect massive binary systems to terminate as spectacular collisions between neutron stars (NSs) and/or black holes (BHs). These events may lead to short GRBs (Eichler et al. 1989) and thereby to chemical enrichment of the interstellar medium by heavy r-process elements (e.g., Rosswog 2015; Just et al. 2015), aside from the powerful emission of gravitational waves (GWs) as recently detected by LIGO (Abbott et al. 2016a, 2016b). For a general review of pre-SN massive star evolution in binaries, see Langer (2012).

Double neutron star (DNS) systems are of special interest since some of them are observable as radio pulsars (Lorimer & Kramer 2004) for several millions, if not billions, of years before

the NSs eventually collide as a result of orbital GW damping, in case their orbital periods are small enough (i.e., less than about 1 day, depending on the orbital eccentricity). Moreover, given that the detected DNS systems represent the endpoint of exotic binary stellar evolution in which progenitor systems have survived multiple stages of mass transfer, one or more common envelope (CE) episodes, as well as two SN explosions (e.g., Flannery & van den Heuvel 1975; Belczynski et al. 2002; Voss & Tauris 2003; Dewi & Pols 2003; Podsiadlowski et al. 2004, 2005; Dewi et al. 2005; Tauris & van den Heuvel 2006), their observed properties are fossil records of their past evolutionary history. Therefore, DNS systems can be used as key probes of binary stellar astrophysics.

Furthermore, the combination of compact objects in relativistic orbits and radio pulsars being ultra-stable clocks allows for unprecedented tests of gravitational theories in the strong-field regime (Damour & Taylor 1992; Stairs 2003; Kramer et al. 2006b; Damour 2009; Kramer & Wex 2009; Freire et al. 2012; Wex 2014; Kaspi & Kramer 2016). Finally, DNS systems help to constrain the long-sought-after equation of state (EoS) of nuclear matter at high densities (Lattimer & Prakash 2007; Kehl et al. 2016; Özel & Freire 2016), although tighter constraints have so far been obtained from massive NSs
with an orbiting white dwarf (WD) companion (Demorest et al. 2010; Antoniadis et al. 2013).

In recent years, the discovery rate of DNS systems has increased, and there are now more than 15 DNS systems known. The coming decade is expected to reveal a large number of discoveries of new DNS systems, as well as their progenitors and merger remnants. The Square Kilometre Array (SKA) is predicted to increase the number of known radio pulsars by a factor of 5 to 10, thus resulting in a total of more than 100 known DNS systems (Keane et al. 2015). The Five-hundred-meter Aperture Spherical Telescope (FAST) is also expected to contribute with a significant number of new radio pulsars (Smits et al. 2009), including new discoveries of DNS binaries. Indirect evidence of DNS mergers is thought to come from detections of short GRBs and their optical afterglows (Berger 2014; Metzger 2017). Moreover, following the first detections of high-frequency GWs from merging double BH binaries (e.g., Abbott et al. 2016a, 2016b), LIGO/VIRGO is also expected to deliver direct evidence of merging DNS systems in the near future (Abadie et al. 2010). Progenitors of DNS systems, the high-mass X-ray binaries (HMXBs), are continuously being discovered with ongoing X-ray missions (INTEGRAL, Swift, XMM, and Chandra; see, e.g., Chaty 2013) and will be supported by upcoming X-ray telescopes like eROSITA, eXTP, and Athena. Hence, we are currently in an epoch where a large wealth of information on NS binaries is available. In light of this, it is important to explore and understand the formation of DNS systems in more detail.

Before presenting the scope of this paper in more detail (Section 1.6), in the following subsections we summarize the background, the current state of the art, and the main motivations for this investigation.

## 1.1. Galactic and Magellanic Cloud Populations of NSs

In our galaxy, the Milky Way, it is estimated that some $10^8–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<fnref n="11" /> (Manchester et al. 2005), there are ~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 ~30 known radio pulsars (Crawford et al. 2001; Ridley et al. 2013), besides an impressive number of ~140 HMXBs (Coe & Kirk 2015; Antoniou & Zezas 2016), some of which are potential progenitors of DNS systems.

## 1.2. Recent DNS Discoveries

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_⊙$), deviating significantly from the “canonical” value of ~1.35 $M_⊙$ (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–1.44 $M_⊙$ 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 (~480 Myr) is on the short side.

## 1.3. Main Characteristics of DNS Systems

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_{NS}$), spin periods ($P$), their time derivatives ($\dot{P}$) and hence surface magnetic fields ($B$), orbital periods ($P_{orb}$), and eccentricities ($e$), as well as kinematic properties such as systemic velocities ($v_{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_{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–185 ms) and their derivatives ($\dot{P}$ = $2.2 \times 10^{-20}–1.8 \times 10^{-17}$ s s⁻¹) are direct measures of the efficiency of the recycling process. The corresponding derived surface dipole $B$-fields are approximately in the range 0.3–18 × 10⁹ G. Table 1 provides the observed ranges of key parameters of DNS systems.

The observed orbital periods now span between $P_{orb}$ = 0.10–45 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_{sys}$ = $v_{LSR}$ ≈ 25–240 km s⁻¹), 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.

<footnote n="11">¹¹ http://www.atnf.csiro.au/research/pulsar/psrcat/</footnote>
<table><caption>**Table 1**
Observed Ranges of Key Properties of DNS Systems</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td colspan="2"><strong>Properties of Recycled (Old) NSs:</strong></td></tr><tr><td>Spin period, <i>P</i></td><td>23–185 ms</td></tr><tr><td>Period derivative, <i>˙P</i></td><td>(0.027–18) × 10<sup>-18</sup> s s<sup>-1</sup></td></tr><tr><td>Surface dipole <i>B</i>-field, <i>B</i></td><td>(0.29–18) × 10<sup>9</sup> G</td></tr><tr><td>Mass, <i>M</i><sub>NS,1</sub></td><td>1.32–1.56 <i>M</i><sub>☉</sub><sup>a</sup></td></tr><tr><td colspan="2"><strong>Properties of Young NSs:</strong></td></tr><tr><td>Spin period, <i>P</i></td><td>144–2773 ms</td></tr><tr><td>Period derivative, <i>˙P</i></td><td>(0.89–20) × 10<sup>-15</sup> s s<sup>-1</sup></td></tr><tr><td>Surface dipole <i>B</i>-field, <i>B</i></td><td>(2.7–5.3) × 10<sup>11</sup> G</td></tr><tr><td>Mass, <i>M</i><sub>NS,2</sub></td><td>1.17–1.39 <i>M</i><sub>☉</sub></td></tr><tr><td colspan="2"><strong>Orbital Properties:</strong></td></tr><tr><td>Orbital period, <i>P</i><sub>orb</sub></td><td>0.10–45 days</td></tr><tr><td>Eccentricity, <i>e</i></td><td>0.085–0.83</td></tr><tr><td>Merger time, <i>τ</i><sub>gwr</sub></td><td>86 Myr → ∞</td></tr><tr><td>Systemic velocity, <i>v</i><sub>sys</sub></td><td>25–240 km s<sup>-1</sup></td></tr></tbody></table>

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_{sys} = v^{LSR}$, is quoted with respect to the local standard of rest (Section 2.2).

<sup>a</sup> 1.32 $M_{☉}$. Mark an upper limit to the lowest mass of the first-born NS.

## 1.4. Résumé of DNS Formation

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<fnref n="12" /> 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<fnref n="13" /> 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

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Illustration of the formation of a DNS system that merges within a Hubble time and produces a single BH, following a powerful burst of GWs and a short GRB. Acronyms used in this figure—ZAMS: zero-age main sequence; RLO: Roche-lobe overflow (mass transfer); He-star: helium star; SN: supernova; NS: neutron star; HMXB: high-mass X-ray binary; CE: common envelope; BH: black hole.</figcaption></figure>

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).

## 1.5. Major Uncertainties in DNS Formation

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

<footnote n="12">¹² See brief discussion given in Section 4.2 for an alternative “double core scenario” (Brown 1995; Dewi et al. 2006) in which CE evolution with a NS is avoided.</footnote>

<footnote n="13">¹³ The secondary (initially less massive) star may be a 5–7 $M_{☉}$ star which accretes mass from the primary (initially more massive) star to reach the threshold limit for core collapse at ~8–12 $M_{☉}$ (Jones et al. 2013; Woosley & Heger 2015; see also Section 3.1).</footnote>
(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

A CE phase (see Taam & Sandquist 2000; Ivanova et al. 2013, for reviews) develops when the donor star in a HMXB fills its Roche lobe. The large mass ratio between the OB-star donor and the accreting NS causes the orbit to shrink significantly upon mass transfer, whereby the NS is captured inside the envelope of the donor star. As a consequence of the resulting drag force acting on the orbiting NS, efficient loss of orbital angular momentum often leads to a huge reduction in orbital separation prior to the second SN explosion, thus explaining the tight orbits of the observed DNS systems.

In a recent study on CE evolution with massive stars, it was demonstrated by Kruckow et al. (2016) that an inspiralling NS may indeed be able to eject the envelope of its massive star companion, at least from an energy budget point of view. However, they also showed that it is difficult to predict the final post-CE separation for several reasons. First of all, it is difficult to locate the bifurcation point (Tauris & Dewi 2001) separating the ejected envelope from the remaining core within the massive star. Second, while additional energy sources, like accretion energy, may play a significant role in helping facilitate the CE-ejection process, models of time-dependent energy transport in the convective envelope are needed to quantify this. Finally, the effect of the inflated envelopes of the remaining helium cores—as a consequence of the stellar luminosity reaching the Eddington limit in their interiors (e.g., Sanyal et al. 2015)—makes it non-trivial to map the core size of a massive donor star model to the size of the naked helium core (Wolf-Rayet star) left behind. Until future 3D hydrodynamical simulations succeed in ejecting the CE of massive stars, the estimated post-CE separation (and thus the LIGO/VIRGO detection rates from the population synthesis of merging BH/NS binaries) will remain highly uncertain.

### 1.5.2. NS Kicks

The second major remaining issue that needs to be solved is the magnitude and direction of the momentum kick added to a NS during the SN explosion (Janka 2012). In particular, for the application to population synthesis modeling of DNS formation, it is important to identify any differences between the first and the second SN explosion (e.g., Bray & Eldridge 2016 and references therein). Whereas there is ample observational evidence for large NS kicks (typically 400–500 km s⁻¹) in observations of young radio pulsars (Gott et al. 1970; Cordes et al. 1993; Lyne & Lorimer 1994; Kaspi et al. 1996; Hobbs et al. 2005; Chatterjee et al. 2005), it also seems evident that the second SN explosion forming DNS systems involves, on average, significantly smaller kicks (Podsiadlowski et al. 2004; van den Heuvel 2004; Schwab et al. 2010; Beniamini & Piran 2016). Furthermore, there is evidence from observations of pulsars in globular clusters (which have small escape velocities of less than 40 km s⁻¹) that some of them are born with very small kicks. This could suggest a bimodal kick velocity distribution, which might be connected to a bimodality in the

formation mechanism of NSs. This picture is supported by observations of HMXBs (Pfahl et al. 2002).

Due to the continuous growth of supercomputing power and the development of efficient and highly parallelized numerical simulation tools, 2D and 3D simulations of stellar core collapse and explosions with sophisticated treatment of the complex microphysics have become feasible and demonstrate the viability of the neutrino-driven and MHD explosion mechanisms in principle (Janka 2012). Initiating the SN blast wave with an approximate neutrino transport description in 2D and 3D simulations, it could be demonstrated (Wongwathanarat et al. 2013) that mass-ejection asymmetries imprinted by hydrodynamic instabilities in the SN core during the initiation of the explosion can lead to NS (and BH) kicks compatible with the measured velocities of young pulsars. Moreover, a first systematic attempt to establish the progenitor-explosion-remnant connection based on the neutrino-driven mechanism has revealed the interesting possibility that the NS-BH transition might occur in stars below $20 M_\odot$ (Ugliano et al. 2012), which seems to be in line with conclusions drawn from observational SN progenitor associations (Smartt 2009). This result should have consequences for the predicted ratio of the number of compact objects (NSs versus BHs) of the GW mergers that LIGO/VIRGO will detect.

It has recently been demonstrated (Tauris et al. 2013, 2015) that close-orbit DNS systems (i.e., those that are tight enough to merge within a Hubble time due to GWs) must form via ultra-stripped SNe when the last star explodes. The reason being that Case BB RLO, i.e., mass transfer via Roche-lobe overflow from a naked helium star in the post-HMXB/post-CE system, causes the NS to significantly strip its evolved helium star companion almost to a naked metal core prior to its explosion. This has an important effect on the number and properties of surviving DNS systems—in particular, in terms of their kinematic properties—as we shall discuss rigorously in more detail.

To learn about the progenitor binaries that formed DNS systems and the conditions of the SN explosions that produced the second-born NSs, many studies in the literature have analyzed or simulated a large number of SNe and compared the outcome with observations of DNS systems (Flannery & van den Heuvel 1975; Hills 1983; Kornilov & Lipunov 1984; Radhakrishnan & Shukre 1985; Dewey & Cordes 1987; Brandt & Podsiadlowski 1995; Kalogera 1996; Tauris & Takens 1998; Wex et al. 2000; Andrews et al. 2015). We therefore conclude our investigation presented in this paper with an extended analysis of the dynamical impact of SNe, taking into account the latest observational constraints on many DNS systems. Our analysis leads to firm results regarding the kick magnitude and direction, as well as the mass of the exploding star in the second SN events.

### 1.5.3. Mass Distribution of NSs

A third main aspect of DNS systems which is far from being understood—and possibly related to the above discussion on SNe and binary stellar evolution in general—is the mass distribution of NSs, and why it differs from the NSs measured in systems with WD companions. Antoniadis et al. (2017) recently analyzed the mass distribution of 32 MSPs in orbit with (mostly) WD companions. They found evidence for a bimodal mass distribution with a low-mass component centered at $1.39 \pm 0.03 M_\odot$ and dispersed by $0.06 M_\odot$, and a high-mass component with a mean
of about $1.81 \pm 0.10 M_\odot$ and dispersed by $0.18 M_\odot$. The diversity in spin and orbital properties of high-mass NSs suggests that this mass difference is most likely not a result of the recycling process, but rather reflects differences in the NS birth masses. For the case of DNS systems, we investigate in this work the total amount of mass accreted by the first-born NS at various phases of evolution.

## 1.6. The Structure and Aims of this Paper

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_{\text{orb}}$, $P$) and ($P_{\text{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.

## 2. Observational Data of DNS Systems

In Table 2, we list 15 DNS systems, including a couple of sources for which the DNS nature is not confirmed. The vast majority of these DNS systems are found in the Galactic disk; only two sources are found in a globular cluster. These two latter sources will be disregarded for further discussions in this paper related to the recycling process and impact of the second SN explosion. The reason is that these two systems almost certainly formed by secondary exchange encounters, i.e., where already-recycled pulsars change stellar companions or are being exchanged into new binaries because of close encounters. In the cases of PSRs B2127+11C and B1807–2500B, such arguments were given by Prince et al. (1991), Lynch et al. (2012), and Verbunt & Freire (2014). In these dynamical processes, information on all traces of the past evolutionary links to their former companions is

lost and cannot be recovered from observations of the present binary systems.

In Table 2, we indicate whether an observed NS is a recycled pulsar or a young (second-born) NS component. Only in the PSR J0737−3039 system are both NS components detected as radio pulsars. This system is therefore known as the *double pulsar*.

### 2.1. Selection Effects

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 \text{ cm}^{-3} \text{ 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_{\text{obs}}$) is more than a few percent of the orbital period ($P_{\text{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_{\text{obs}}/P_{\text{orb}} \ge 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.
<table><caption>**Table 2**

Properties of 15 DNS Systems with Published Data (Including a Few Unconfirmed Candidates)</caption><thead><tr><th>Radio Pulsar</th><th>Type</th><th>P (ms)</th><th>P&#775; (10<sup>-18</sup>)</th><th>B (10<sup>9</sup> G)</th><th>P<sub>orb</sub> (days)</th><th>e</th><th>M<sub>psr</sub> (M<sub>☉</sub>)</th><th>M<sub>comp</sub> (M<sub>☉</sub>)</th><th>M<sub>total</sub> (M<sub>☉</sub>)</th><th>&delta; (deg)</th><th>Dist. (kpc)</th><th>v<sup>LSR</sup> (km s<sup>-1</sup>)</th><th>&tau;<sub>gwr</sub> (Myr)</th></tr></thead><tbody><tr><td>J0453+1559 (1)</td><td>recycled</td><td>45.8</td><td>0.186</td><td>0.92</td><td>4.072</td><td>0.113</td><td>1.559</td><td>1.174</td><td>2.734</td><td>&mdash;</td><td>1.07</td><td>52</td><td>&infin;</td></tr><tr><td>J0737-3039A (2)</td><td>recycled</td><td>22.7</td><td>1.76</td><td>2.0</td><td>0.102</td><td>0.088</td><td>1.338</td><td>1.249</td><td>2.587</td><td>&lt;3.2</td><td>1.15</td><td>33</td><td>86</td></tr><tr><td>J0737-3039B (2)</td><td>young</td><td>2773.5</td><td>892</td><td>490</td><td>&mdash;||&mdash;</td><td>&mdash;||&mdash;</td><td>1.249</td><td>1.338</td><td>&mdash;||&mdash;</td><td>130 &plusmn; 1</td><td>&mdash;||&mdash;</td><td>&mdash;||&mdash;</td><td>&mdash;||&mdash;</td></tr><tr><td>J1518+4904 (3)</td><td>recycled</td><td>40.9</td><td>0.0272</td><td>0.29</td><td>8.634</td><td>0.249</td><td>c</td><td>c</td><td>2.718</td><td>&mdash;</td><td>0.63</td><td>28</td><td>&infin;</td></tr><tr><td>B1534+12 (4)</td><td>recycled</td><td>37.9</td><td>2.42</td><td>3.0</td><td>0.421</td><td>0.274</td><td>1.333</td><td>1.346</td><td>2.678</td><td>27 &plusmn; 3</td><td>1.05</td><td>136</td><td>2730</td></tr><tr><td>J1753-2240 (5)</td><td>recycled</td><td>95.1</td><td>0.970</td><td>2.7</td><td>13.638</td><td>0.304</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>3.46</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>J1755-2550 (6)<sup>a</sup></td><td>young</td><td>315.2</td><td>2430</td><td>270</td><td>9.696</td><td>0.089</td><td>&mdash;</td><td>&gt;0.40</td><td>&mdash;</td><td>&mdash;</td><td>10.3</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>J1756-2251 (7)</td><td>recycled</td><td>28.5</td><td>1.02</td><td>1.7</td><td>0.320</td><td>0.181</td><td>1.341</td><td>1.230</td><td>2.570</td><td>&lt;34</td><td>0.73</td><td>41<sup>d</sup></td><td>1660</td></tr><tr><td>J1811-1736 (8)</td><td>recycled</td><td>104.2</td><td>0.901</td><td>3.0</td><td>18.779</td><td>0.828</td><td>&lt;1.64</td><td>&gt;0.93</td><td>2.57</td><td>&mdash;</td><td>5.93</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>J1829+2456 (9)</td><td>recycled</td><td>41.0</td><td>0.0525</td><td>0.46</td><td>1.176</td><td>0.139</td><td>&lt;1.38</td><td>&gt;1.22</td><td>2.59</td><td>&mdash;</td><td>0.74</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>J1906+0746 (10)<sup>a</sup></td><td>young</td><td>144.1</td><td>20300</td><td>530</td><td>0.166</td><td>0.085</td><td>1.291</td><td>1.322</td><td>2.613</td><td>&mdash;</td><td>7.40</td><td>&mdash;</td><td>309</td></tr><tr><td>J1913+1102 (11)</td><td>recycled</td><td>27.3</td><td>0.161</td><td>0.63</td><td>0.206</td><td>0.090</td><td>&lt;1.84</td><td>&gt;1.04</td><td>2.88</td><td>&mdash;</td><td>&mdash;</td><td>&mdash;</td><td>~480</td></tr><tr><td>B1913+16 (12)</td><td>recycled</td><td>59.0</td><td>8.63</td><td>7.0</td><td>0.323</td><td>0.617</td><td>1.440</td><td>1.389</td><td>2.828</td><td>18 &plusmn; 6</td><td>9.80</td><td>240</td><td>301</td></tr><tr><td>J1930-1852 (13)</td><td>recycled</td><td>185.5</td><td>18.0</td><td>18</td><td>45.060</td><td>0.399</td><td>&lt;1.32</td><td>&gt;1.30</td><td>2.59</td><td>&mdash;</td><td>1.5</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>J1807-2500B (14)<sup>a</sup></td><td>GC</td><td>4.2</td><td>0.0823</td><td>0.18</td><td>9.957</td><td>0.747</td><td>1.366</td><td>1.206</td><td>2.572</td><td>&mdash;</td><td>3.0</td><td>&mdash;</td><td>&infin;</td></tr><tr><td>B2127+11C (15)</td><td>GC</td><td>30.5</td><td>4.99</td><td>3.8</td><td>0.335</td><td>0.681</td><td>1.358</td><td>1.354</td><td>2.713</td><td>&mdash;</td><td>12.9</td><td>&mdash;</td><td>217</td></tr></tbody></table>

**Notes.** *Ṗ* values are not corrected for the (in most cases negligible) Shklovskii effect (Shklovskii 1970). *B*-field values are estimated from Equation (1). All the NS masses quoted with four significant figures have uncertainties of ≲0.010 *M*<sub>☉</sub>. All values of τ<sub>gwr</sub> > 50 Gyr are shown with ∞.

<sup>a</sup> Not a confirmed DNS system—could also be a WD+NS binary.

<sup>b</sup> Based on the median value of *v*<sup>LSR</sup>; see Section 2.2 and Table 3.

<sup>c</sup> See Section 8.3 for updated NS masses of PSR J1518+4904.

<sup>d</sup> Based on the 2D proper motions estimated in Ferdman et al. (2014).

**References.** (1) Martinez et al. (2015). (2) Kramer et al. (2006b), Breton et al. (2008), Ferdman et al. (2013). (3) Janssen et al. (2008). (4) Fonseca et al. (2014). (5) Keith et al. (2009). (6) Ng et al. (2015), C. Ng et al. (2017, in preparation). (7) Faulkner et al. (2005), Ferdman et al. (2014). (8) Corongiu et al. (2007). (9) Champion et al. (2004). (10) Lorimer et al. (2006), van Leeuwen et al. (2015). (11) Lazarus et al. (2016). (12) Hulse & Taylor (1975), Kramer (1998), Weisberg et al. (2010). (13) Swiggum et al. (2015). (14) Lynch et al. (2012). (15) Anderson et al. (1990); Jacoby et al. (2006).

(the shorter the integration, the smaller is the bias) and the particular search algorithms being used in each survey.

## 2.2. Estimating the Systemic Velocities of DNS Systems Relative to Their LSR

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*<sub>*T*</sub><sup>SSB</sup>, 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*<sub>*R*</sub><sup>SSB</sup> (based on the calculated *v*<sub>*T*</sub><sup>SSB</sup>): *v*<sub>*R*</sub><sup>SSB</sup> =

*v*<sub>*T*</sub><sup>SSB</sup> cot χ, 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*<sup>LSR</sup>, 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σ and 2σ ranges, as well as the 90% confidence limit for the upper limit of the 3D systemic velocity, *v*<sup>LSR,max</sup>, based on our MC simulations (strictly speaking, since the *v*<sup>LSR</sup> distributions are non-Gaussian, these 1σ and 2σ 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<sup>14</sup> in our analysis of DNS kinematics discussed further in this work.

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

<footnote n="14"><sup>14</sup> In the case of PSR B1913+16, we adopted a default distance based on H I absorption measurements (Weisberg et al. 2008).</footnote>
<table><caption>**Table 3**
Estimated 3D Systemic LSR Velocities for Galactic Disk DNS Systems with Proper Motion Measurements</caption><thead><tr><th>Radio Pulsar</th><th><i>l</i> (deg)</th><th><i>b</i> (deg)</th><th>Distance (kpc)</th><th>μ<sub>α</sub> (mas yr<sup>-1</sup>)</th><th>μ<sub>δ</sub> (mas yr<sup>-1</sup>)</th><th><i>v</i><sub>T</sub><sup>SSB</sup> (km s<sup>-1</sup>)</th><th><i>v</i><sub>T</sub><sup>LSR</sup> (km s<sup>-1</sup>)</th><th><i>v</i><sub>1σ</sub><sup>LSR</sup> (km s<sup>-1</sup>)</th><th><i>v</i><sub>2σ</sub><sup>LSR</sup> (km s<sup>-1</sup>)</th><th><i>v</i><sub>90% CL</sub><sup>LSR, max</sup> (km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>J0453+1559</td><td>184.12</td><td>-17.14</td><td>1.07</td><td>-5.5</td><td>-6.0</td><td>41 ± 8</td><td>52</td><td>36-85</td><td>26-197</td><td>102</td></tr><tr><td>J0737-3039A</td><td>245.24</td><td>-4.50</td><td>1.15 (a)</td><td>-3.8</td><td>2.1</td><td>25 ± 4</td><td>33</td><td>18-57</td><td>10-124</td><td>68</td></tr><tr><td>J1518+4904</td><td>80.81</td><td>54.28</td><td>0.63</td><td>-0.7</td><td>-8.5</td><td>26 ± 5</td><td>28</td><td>20-49</td><td>15-123</td><td>61</td></tr><tr><td>B1534+12</td><td>19.85</td><td>48.34</td><td>1.05 (b)</td><td>1.5</td><td>-25.3</td><td>127 ± 1</td><td>136</td><td>118-224</td><td>115-512</td><td>276</td></tr><tr><td>J1756-2251</td><td>6.50</td><td>0.95</td><td>0.73 (c)</td><td>-2.4</td><td>5.5</td><td>21 ± 7</td><td>41</td><td>21-71</td><td>11-123</td><td>82</td></tr><tr><td>B1913+16</td><td>49.97</td><td>2.12</td><td>9.80 (d)</td><td>-1.2</td><td>-0.8</td><td>69 ± 22</td><td>240</td><td>174-273</td><td>83-401</td><td>290</td></tr><tr><td>J0453+1559</td><td>184.12</td><td>-17.14</td><td>0.52</td><td>-5.5</td><td>-6.0</td><td>20 ± 4</td><td>27</td><td>20-43</td><td>15-99</td><td>52</td></tr><tr><td>J1518+4904</td><td>80.81</td><td>54.28</td><td>0.96</td><td>-0.7</td><td>-8.5</td><td>39 ± 8</td><td>39</td><td>28-70</td><td>21-184</td><td>89</td></tr><tr><td>B1913+16</td><td>49.97</td><td>2.12</td><td>5.25</td><td>-1.2</td><td>-0.8</td><td>37 ± 7</td><td>132</td><td>98-168</td><td>67-229</td><td>179</td></tr></tbody></table>

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.

## 2.3. Derived Surface B-fields

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 \approx 3.2 \times 10^{19} 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 $j \times 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 $\alpha$, Tauris et al. (2012) obtained$^{15}$

$$ 
\begin{aligned}
B &= \sqrt{\frac{c^3 I P \dot{P}}{4\pi^2 R^6} \frac{1}{1 + \sin^2 \alpha}} \\
&\approx 1.3 \times 10^{19} G \sqrt{P\dot{P}} \left(\frac{M}{1.4 M_{\odot}}\right)^{3/2} \sqrt{\frac{1}{1 + \sin^2 \alpha}}.
\end{aligned} 
\quad (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 \approx 1.0 \times 10^{19} 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.

## 2.4. Measuring NS Masses

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

<figure><img src="image_4.png" /><figcaption>**Figure 2.** ($P, \dot{P}$) diagram of all radio pulsars detected in DNS systems. Data are given in Table 2. The solid gray lines represent constant surface dipole B-fields and the dotted black lines represent constant characteristic ages. The $\dot{P}$ values are not corrected for the (mostly negligible) Shklovskii effect (Shklovskii 1970). For the globular cluster pulsars, $\dot{P}$ must also be corrected for the (potentially significant) NS acceleration in the cluster potential.</figcaption></figure>

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^{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}_{\text{orb}}$).

To further check our results against the sensitivity of the DM distance, we recalculated the values of $v^{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}^{LSR}$ and $v_{2\sigma}^{LSR}$ only changed by a few km s$^{-1}$ for PSRs J0453+1559 and J1518+4904, whereas for PSR B1913+16 these interval ranges widened on each end by ~18 km s$^{-1}$ and ~12 km 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.

<footnote n="15">15 Using slightly larger values for the NS radius and moment of inertia than the standard values of $R = 10$ km and $I = 10^{45}$ g cm$^2$.</footnote>
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,
ω, 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 (ω, γ) 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, P̂<sub>orb</sub>.
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*<sub>3</sub> and *ζ* in the Freire & Wex (2010)
parameterization. These allow for a direct determination of the
mass of the pulsar's companion (r = GM<sub>comp</sub>/c<sup>2</sup>) 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 P̂<sub>orb</sub> 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 ω.

## 2.5. Observed Distribution of Masses in DNS Systems

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 <~0.01 M<sub>☉</sub>, whereas the error bars for the pulsar masses in the NS+WD systems plotted here range between 0.1%–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 *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 ~0.1 M<sub>☉</sub>. 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.

## 3. Evolution of OB-star Binaries from the ZAMS to the HMXB Stage

In the following section, we discuss the early stages of
massive binary star evolution and, in particular, the observed
Galactic population of HMXBs. We examine their destiny with
the aim of understanding the evolutionary links between
HMXBs and DNS systems.

3.1. Pre-HMXB Evolution

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
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Distribution of NS star masses in DNS systems (top, central) and NS +WD systems (bottom panel). Recycled pulsars are indicated by (r) and blue or yellow color; young pulsars are indicated by (y) and green or red color. It is evident that NSs that formed second in these binaries (i.e., the non-recycled, young pulsars) in general have a smaller mass than the recycled NSs. See Section 2.5. Note that none of the colored bars overlap—all NS data are visible.</figcaption></figure>

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 (≈10 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–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.

## 3.2. Observations of HMXB Systems

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$ 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.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Extended Corbet diagram showing the different populations of HMXBs (X-ray pulsars) and DNS systems (recycled radio pulsars) with measured values of both $P_{\text{orb}}$ and P. A linear correlation among the DNS systems is given by Equation (6). Interesting questions are which, and how, HMXBs evolve to become DNS systems—see the text for discussions, and Section 3.3 for information on the data.</figcaption></figure>

### 3.2.2. sgHMXBs

sgHMXBs host a supergiant star of spectral type O8–B1 I/II, characterized by an intense, dense, and highly supersonic radiatively-driven stellar wind (e.g., Martínez-Núñez et al. 2017). B-type supergiants typically possess slower winds ($v_\infty = 1000 \text{ km s}^{-1}$) than O-type supergiants ($v_\infty = 2000 \text{ km s}^{-1}$; Negueruela 2010). There are $\sim 16$ so-called “classical” persistent sgHMXBs, most of them being close systems with a compact object in a short and circular orbit, directly accreting from the stellar wind through a Bondi–Hoyle-like accretion process (see Section 4.1). Such wind-fed systems exhibit a luminous and persistent X-ray emission ($L_X = 10^{36-38} \text{ erg s}^{-1}$), with superimposed large variations on short timescales and a cutoff (10–30 keV) power-law X-ray spectrum. Located in the upper-left part of the Corbet diagram (small orbital periods of $P_{\text{orb}} \simeq 3–10$ days and long spin periods of $P \simeq 100–10,000$ s; Figure 4), sgHMXBs do not show any correlation between $P_{\text{orb}}$ and P, since there is no net transfer of angular momentum via wind accretion (Waters & van Kerkwijk 1989). The detection of regular pulsations implies that they host young NSs with surface magnetic fields of $B \sim 10^{11–12}$ G. Nearly half of the sgHMXBs exhibit a substantial intrinsic, local extinction $n_H \gtrsim 10^{22} \text{ cm}^{-2}$, with a compact object deeply embedded in the dense stellar wind (such as the highly obscured IGR J16318-4848; Chaty & Rahoui 2012). Likely in transition to RLO-HMXBs, these systems are characterized by slow winds which, depending on the degree of flow through the inner Lagrangian point (L1), may cause shrinking of the orbit and the spiral-in of the compact object over time, leading to a CE phase (Sections 1.5.1 and 4.2). An example of a sgHMXB where accretion takes place both through RLO and stellar wind is the BH system Cyg X-1.

A significant subclass of sgHMXBs comprises some $\sim 12$ (and a few candidate) supergiant fast X-ray transients (SFXTs; see Negueruela et al. 2006; Walter et al. 2015). These systems, characterized by a compact object orbiting with $P_{\text{orb}} \simeq 3–100$ days in a circular or eccentric orbit, and typically with $P \simeq 100–1000$ s, span a large area in the Corbet diagram,
mostly in between BeHMXBs and sgHMXBs (Figure 4). They
are not persistent sources, but exhibit short and intense X-ray
outbursts, an unusual characteristic among HMXBs, rising in
tens of minutes up to a peak luminosity $L_X \simeq 10^{35-37}$ erg s$^{-1}$,
lasting for a few hours, and alternating with long (~70 days)
epochs of quiescence at $L_X \sim 10^{32-34}$ erg s$^{-1}$, with an
impressive variability factor $L_{max}/L_{min}$ going up to $10^2-10^6$
(Romano et al. 2015). Various processes have been invoked to
explain these flares, such as wind inhomogeneities, magneto/
centrifugal accretion barriers, transitory accretion disks, etc.
(see, e.g., Chaty 2013; Walter et al. 2015, and references
therein). To explain the distribution of sgHMXBs in the
Corbet diagram, Liu et al. (2011) argued that the majority of
the systems evolved directly from OB-type main-sequence
star-NS systems (i.e., without a significant accretion history
since the formation of the NS), whereas a minority of the
systems evolved from BeHMXBs (i.e., SFXTs with an
accretion history and thus located within the area of
BeHMXBs).

3.2.3. RLO-HMXBs

Beginning (atmospheric) RLO-HMXBs host a massive star
filling its Roche lobe, where accreted matter flows via the inner
Lagrangian point (L1) to form an accretion disk (similarly to
the situation in LMXBs). This process is often referred to as
beginning or “atmospheric RLO” (Savonije 1978, 1979, 1983).
They constitute the classical bright pulsing HMXBs (i.e.,
Cen X-3, SMC X-1, and LMC X-4), which have high X-ray
luminosities ($L_X \sim 10^{38}$ erg s$^{-1}$) during outbursts (possibly
caused by accretion disk instabilities; Dubus et al. 2001). There
are only a few such sources, located in the lower-left part of the
Corbet diagram (Figure 4), characterized by short orbital and
spin periods. The lifetime of these systems before they become
dynamically unstable is expected to be relatively short (at most
a few $10^4$ yr; Savonije 1978).

3.3. HMXB Population Statistics

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 ~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 kpc of an SFC, with an
intercluster distance of 1.7 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 ~$10^7$ 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.

3.4. Evolution and Destiny of HMXB Systems

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<fnref n="16" /> (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*

If the orbital period of the binary at the time of RLO is
relatively short ($\lesssim 1$ yr; Taam et al. 1978; Terman et al. 1995;
Taam 1996), then the binding energy of the donor star envelope
remains too large to allow for a successful ejection via the
orbital energy released in the spiral-in (e.g., Kruckow et al.
2016). In this case, the NS is expected to spiral to the center,
leading to the complete coalescence of the system. The
resulting objects are referred to as Thorne–Żytkow objects
(TŻOs; Thorne & Żytkow 1975, 1977). These objects will

<footnote n="16">¹⁶ A NS orbiting a massive star induces a tide in the massive star, which will try to spin it up and ultimately make it co-rotate with the NS’s orbit. However, for a sufficiently large mass ratio, there is not enough angular momentum in the NS’s orbit to bring the massive star into co-rotation. The transfer of angular momentum from the orbit to the star causes continued shrinking of the orbit.</footnote>
appear as very cool red supergiants. Since these are difficult to distinguish from normal red supergiants, it is presently not clear whether they actually exist. Since this is the possible fate for the far majority of known HMXBs shown in Figure 4, their birth rate in the Galaxy is expected to be quite high (~2 × 10⁻⁴ yr⁻¹; Podsiadlowski et al. 1995). Depending on the uncertain lifetime of this phase (which is limited, e.g., by wind mass loss), a few to 10% of all red supergiants with a luminosity comparable to or above the Eddington limit for a NS ($L_{Edd} \sim 10^5 L_{\odot}$) may harbor NS cores. Massive TŻOs would be distinguishable from normal red supergiants through their anomalously large abundances of proton-rich elements, in particular molybdenum (Biehle 1991; Cannon 1993). Despite numerous searches, to date only one candidate TŻO has been identified (HV 2112 in the SMC; Levesque et al. 2014) and its interpretation as a TŻO has remained controversial (macfarlane & de Mink 2016). This suggests that the lifetime of the TŻO phase (if it exists) is much shorter than previously estimated.

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

If the orbital period of the HMXB hosting a NS is relatively long (≥1 yr), the binding energy of the donor star envelope at the onset of the CE will be sufficiently small to allow envelope ejection, depending on the hydrodynamical conversion of the released orbital energy and the mass of the donor star. Kruckow et al. (2016) demonstrated that solutions exist for ejecting the CE in such post-HMXB systems with donor star masses below ~22 $M_{\odot}$ and ~26 $M_{\odot}$, for metallicities equivalent to that of the Milky Way average and $Z = Z_{\odot}/50$, respectively. From the orbital period distribution of the known Galactic HMXBs, we conclude that very few, if any, of these binaries will produce DNS systems, a fact that is not in contradiction with population statistics given that the radio lifetime of DNS systems ($10^8-10^{10}$ yr for the recycled pulsar; see Figure 2) is much longer than the X-ray lifetime of HMXBs (typically between some $10^5$ and a few $10^6$ yr). For recycled pulsars, the characteristic age is an approximate measure of the remaining lifetime as an active radio pulsar, given that $\tau = P/(2\dot{P}) = E_{\text{rot}}/|\dot{E}_{\text{rot}}|$ (for discussions, see Tauris et al. 2012). The requirement of a wide orbit for the HMXB system to survive the CE phase has interesting consequences for the magnitude of the NS kick in the first SN explosion; see discussion in Section 6.9.

In the case of successful CE ejection, the post-CE system will be a much closer binary consisting of the NS in orbit with the hydrogen-exhausted core of the massive secondary star. If this remaining helium star has mass $≥3-4 M_{\odot}$, it will have significant stellar wind (e.g., Yoon et al. 2010 and references therein). If even a small fraction of this wind is accreted by the NS (see Section 4.3), the system will again appear as an X-ray binary. Eventually, and in many cases following an additional phase of mass transfer (Case BB RLO, Section 4.4), the helium star will explode in an (ultra-)stripped SN of Type Ib/Ic and will itself produce a NS (Figure 1). Depending on the magnitude of the natal NS kick, the system may be disrupted in this second SN event (see Sections 6.5 and 6.9). If this is the case, both NSs (one a young pulsar, the other a relatively old and mildly recycled NS) will move apart as runaway NSs (Radhakrishnan & Shukre 1985; Tauris & Takens 1998; Lorimer et al. 2004). On the other hand, if the system remains bound, the surviving binary finally becomes a DNS system.

### 3.5. ULXs: Origin and Connection to DNS Progenitors?

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–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_2$ 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)

$$ \frac{a}{a_0} = \left( \frac{q_0(1 - \beta) + 1}{q(1 - \beta) + 1} \right)^{\frac{3\beta - 5}{1 - \beta}} \left( \frac{q_0 + 1}{q + 1} \right) \left( \frac{q_0}{q} \right)^2, \quad (2) $$

where $a_0$ and $a$ refer to the orbital separation before and during RLO, respectively, and where $q_0$ and $q$ represent the mass ratios at these two epochs. The parameter $\beta$ 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.

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 (~1.8 × $10^{-8} M_{\odot}$ yr⁻¹), 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
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Orbital evolution of IMXBs and HMXBs based on the isotropic re-emission model. The plot shows the decrease in orbital separation (in units of the pre-RLO orbital separation, $a_0$) as a function of donor star mass, for an IMXB system ($M_2 = 6 M_\odot$) and a HMXB system ($M_2 = 12 M_\odot$). Minimum values of ($a/a_0$) are shown with circles for different values of $\beta$ in steps of 0.2. NS ULXs are likely to be IMXBs where the RLO can remain stable for a relatively long time, unlike the situation for HMXBs which evolve with strong orbital decay—see Section 3.5 for discussions.</figcaption></figure>

above the Eddington limit (even on the order of ~100 $\dot{M}_{\text{Edd}}$; Israel et al. 2017a). However, given that IMXBs and HMXBs immediately reach values of $|\dot{M}_2| > 10^{-5} M_\odot \text{yr}^{-1}$ after the onset of RLO, the estimate of $\beta > 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 > 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-10^4 \text{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| > 10^{-5} M_\odot \text{yr}^{-1}$ for a few $10^5-10^6 \text{yr}$ (e.g., Tauris et al. 2000, 2011), followed by RLO at a lower rate for some $10^6-10^7 \text{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.

We have recently been notified that optical monitoring of the NS ULX system NGC 7793 P13 (Motch et al. 2014) has revealed a ~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

<table><caption>**Table 4**
Linking IMXBs, HMXBs, NS ULXs, and DNS Systems</caption><thead><tr><th>M<sub>2</sub>/M<sub>☉</sub></th><th>Stability</th><th>NS ULX source</th><th>Final outcome</th></tr></thead><tbody><tr><td>≤5</td><td>stable RLO</td><td>yes</td><td>NS + He/CO WD</td></tr><tr><td>6–9</td><td>DDI → CE</td><td>?</td><td>NS + CO/ONeMg WD</td></tr><tr><td>≥10</td><td>CE</td><td>no</td><td>DNS</td></tr></tbody></table>

**Note.** The division between IMXBs and HMXBs is somewhat unclear but typically at a donor star mass of $M_2 \simeq 8-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).

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.

## 4. Amount of Mass Accreted by Recycled Pulsars in DNS Systems

A first estimate of the amount of matter accreted by recycled pulsars in different types of binary pulsar systems was made by Taam & van den Heuvel (1986). Here, to infer the birth masses of the first-born NSs in DNS systems (and allow for direct comparison with the masses of the second-born NSs, or NSs in other binaries like HMXBs), we have identified five consecutive phases of evolution where the first-born NS in a DNS system will potentially accrete material (see Figure 1):

(i) wind accretion in the HMXB stage,

(ii) CE and spiral-in evolution,

(iii) wind accretion from the helium/Wolf-Rayet star,

(iv) Case BB Roche-lobe overflow, and

(v) shell impact from the SN of the secondary star.

We now discuss the expected amount of matter accreted by the NS in each of these five phases.

### 4.1. Wind Accretion in the HMXB Stage

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,wind}| \simeq 10^{-6}-10^{-5} M_\odot \text{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_{wind} \simeq v_{esc} = \sqrt{2GM_2/R_2}$, where $M_2$ and $R_2$ 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_{acc} = 2GM_{NS}/(v_{rel}^2 + c_s^2)$, where $c_s = \sqrt{\gamma P/\rho} \approx 11 \sqrt{T/(10^4 K)}$ km s$^{-1}$ is the local sound speed of the ambient medium ($\gamma$ is the adiabatic index, $P$ is the pressure, $\rho$ is the mass density, and $T$ is the temperature). The relative velocity between the wind and the NS ($v_{rel} = v_{wind} + v_{orb}$) is usually dominated by the wind velocity compared to the orbital velocity, i.e., $v_{wind} \gg v_{orb}$ (except for very tight and/or highly eccentric binaries). Similarly, we have $v_{wind} \gg c_s$. The resulting accretion rate onto the NS, $\dot{M}_{NS} \approx \pi R_{acc}^2 \rho v_{wind}$, can be combined with the continuity equation $\rho = |\dot{M}_{2,wind}|/(4\pi a^2 v_{wind})$ to yield

$$ \dot{M}_{NS} \approx \frac{(GM_{NS})^2 |\dot{M}_{2,wind}|}{a^2 v_{wind}^4}, \quad (3) $$

where $a$ is the orbital separation. An estimate from this equation indicates that, in wind-driven scenarios, $\dot{M}_{NS}$ is typically $10^3-10^4$ times lower than $|\dot{M}_{2,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}_{NS} \approx 10^{-9} M_{\odot} \text{yr}^{-1}$, corresponding to an X-ray luminosity of the order $L_x \approx 10^{37} \text{ erg 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: $\Delta M_{NS} = \dot{M}_{NS} \cdot \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_x \ge 10^{37} \text{ erg 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,wind}|$, is too large. During most of the (main-sequence) lifetime of the HMXB donor star, the value of $|\dot{M}_{2,wind}|$ is often lower than $10^{-6} M_{\odot} \text{yr}^{-1}$. Indeed, from a study of 10 eclipsing HMXBs, Falanga et al. (2015) list typical values of $|\dot{M}_{2,wind}| \approx 10^{-7}-10^{-6} M_{\odot} \text{yr}^{-1}$, based on observations of their average NS accretion rates. Hence, by adopting a value of $\Delta t_2 \sim 10^7 \text{ yr}$ and slightly smaller values of $|\dot{M}_{2,wind}|$ and $\dot{M}_{NS}$, we finally obtain $\Delta M_{NS} \approx \text{a 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).

## 4.2. CE and Spiral-in Evolution

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 $\Delta M_{NS} < 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_{NS} < 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 $\Delta M_{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.

## 4.3. Wind Accretion from the Helium/Wolf–Rayet Star

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., $\Delta M_{NS} \approx |\dot{M}_{He,wind}| \cdot f_{acc} \cdot \Delta t_{He}$, where $\Delta t_{He}$ is the lifetime of the helium star, $|\dot{M}_{He,wind}|$ is its mass-loss rate, and $f_{acc} = \dot{M}_{NS}/|\dot{M}_{He,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_{acc} \simeq \text{a few} \times 10^{-4}$ (see
also Israel 1996), and we obtain values of $\Delta M_{NS} < 4 \times 10^{-4} M_{\odot}$ when integrating throughout the wind accretion phase of the NS–helium star binaries.

### 4.4. Case BB Roche-lobe Overflow

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 ($<10^5$ yr) combined with an Eddington-limited accretion rate onto the NS, which is roughly given by (cgs units in the central term)

$$ \dot{M}_{\text{Edd}} = \frac{4\pi c R_{\text{NS}}}{0.2(1 + X)} \approx 3.0 \times 10^{-8} M_{\odot} \text{ yr}^{-1} \quad (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 $\Delta M_{NS} = 5 \times 10^{-5} - 3 \times 10^{-3} M_{\odot}$, where $\Delta M_{NS}$ is systematically decreasing with increasing values of $P_{\text{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}_{\text{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}_{\text{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-9) \times 10^{-3} M_{\odot}$.

### 4.5. Shell Impact from the SN of the Secondary Star

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 $\Delta M_{NS} \ll 10^{-3} M_{\odot}$.

### 4.6. The NS Accretion Efficiency

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$ ms, in orbit with a low-mass He WD and $P_{\text{orb}} = 10.9$ 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 $\Delta M_{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 $\Delta M_{NS}$ at various stages discussed above may represent overestimates due to magnetospheric and accretion disk effects, which were not taken into account.

### 4.7. Total Amount of Mass Gained by the First-born NS

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 $\Delta M_{NS} < 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.

### 4.8. Spin Periods of Recycled NSs in DNS Systems

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)<fnref n="17" />

$$ \Delta M_{\text{eq}} \approx 0.22 M_{\odot} \frac{(M_{\text{NS}}/M_{\odot})^{1/3}}{P_{\text{ms}}^{4/3}} \quad (5) $$

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

<footnote n="17">¹⁷ Under the assumption that the location of the NS magnetosphere is approximately kept fixed and the accreting NS spins near/at equilibrium while accumulating most of its material.</footnote>
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 ~11 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<sub>☉</sub>) 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}_{NS} \gg \dot{M}_{Edd}$. Such a high accretion rate is apparently evident in an 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-2 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.

## 4.9. Birth Masses of NSs in DNS Systems

The above analysis demonstrates the very limited accretion (<0.02 M<sub>☉</sub>) 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 ~0.1 M<sub>☉</sub>) 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 ($\ge 1.70 M_{☉}$) 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 ~2.0 M<sub>☉</sub>, 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.

## 5. The ($P_{orb}$, $P$, $e$) Correlations

An important diagnostic tool for probing the formation of DNS systems is the investigation of correlations between orbital parameters (orbital period and eccentricity) and the spin period of the recycled pulsars. Any correlation will yield valuable information about the previous mass-transfer phase and the SN explosion. The first such investigations were carried out by Dewi et al. (2002, 2005), Dewi & Pols (2003), Faulkner et al. (2005). A similar method is applied when investigating the formation and evolution of BeHMXBs (Sections 3.2 and 3.4) following the first SN; see, e.g., Pfahl et al. (2002) and Knigge et al. (2011).

In the following, we first demonstrate an empirical correlation between $P_{orb}$ and $P$ (Section 5.1) followed by the theoretical modeling of such a correlation (Section 5.2). We then discuss correlations between $P_{orb}$ and eccentricity (Section 5.3), and $P$ and eccentricity (Section 5.4), with a particular focus on the NS kick from the second SN. More detailed analyses and discussions on NS kicks follow thereafter in Section 6.

### 5.1. The ($P_{orb}$, $P$) Correlation

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

$$P = 44 \text{ ms } (P_{\text{orb}}/\text{days})^{0.26}, \qquad (6)$$

based on a linear regression in the (log $P_{orb}$, 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_{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_{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
<figure><img src="image_8.png" /><figcaption>Figure 6. Spin period as a function of orbital period for all first-born (recycled) NSs in DNS systems. The observational data are plotted with red points (the error bars are too small to be seen). The gray line (Equation (6)) is a fit to the raw data. The green line (Equation (7)) is our estimated empirical relation at birth (i.e., right after the second SN), and the shaded region is its uncertainty.</figcaption></figure>

correlation (green line):

$$P \approx 36 \pm 14 \text{ ms } (P_{\text{orb}}/\text{days})^{0.40 \pm 0.10}. \quad (7)$$

A corresponding theoretical ($P_{\text{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 > 100$ 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$ ms, $\dot{P} = 1.8 \times 10^{-17}$ s s$^{-1}$, $P_{\text{orb}} = 45$ days, and eccentricity $e = 0.40$; see Figure 2). Its discovery is important for testing formation scenarios of such wide-orbit marginally recycled pulsars.

## 5.2. Calculating a Theoretical ($P_{\text{orb}}$, $P$) Correlation

Based on a few assumptions, one can derive a theoretical ($P_{\text{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, $\Delta M_{\text{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

<figure><img src="image_9.png" /><figcaption>Figure 7. Kippenhahn diagram of a 3.0 $M_{\odot}$ helium star undergoing Case BB RLO to a NS in a wide binary with $P_{\text{orb},i} = 50$ days (see data in Table 5). The plot shows the evolving cross-section of the star, in mass coordinates on the y-axis, as a function of the remaining calculated lifetime on the x-axis. The green hatched areas denote zones with convection and the intensity of the blue/ purple color indicates the net energy-production rate. The total mass of the star is shown by the solid black line. For this wide-orbit binary, the onset of the RLO (marked by a red arrow) will occur at a late evolutionary stage (carbon-shell burning) such that the star will undergo core collapse within ~1000 yr, i.e., before the end of the RLO. Mass loss prior to RLO is due to a stellar wind.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 8. Theoretical correlations between the orbital period and spin period of the first-born NSs in DNS systems (blue lines), calculated for three different values of the Eddington accretion efficiency parameter, $X_{\text{Edd}}$, and assuming $M_{\text{He,i}} = 3.0 M_{\odot}$ and no kick during the SN. The observational data are plotted with red points, and the green line is our birth relation based on this data (see Figure 6).</figcaption></figure>

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_{\text{orb}}$, and the eccentricity (see Section 6.3).

Figure 8 shows examples of such theoretical ($P_{\text{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_{\text{Edd}}$. This parameter is a measure of the
<table><caption>**Table 5**
Properties of Case BB RLO Systems with $M_{\text{He,i}} = 3.0 M_{\odot}$</caption><thead><tr><th><i>M</i><sub>He,i</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>orb,i</sub><br>(days)</th><th><i>M</i><sub>He,f</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>core,f</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>orb,f</sub><br>(days)</th><th>Δ<i>M</i><sub>NS</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>NS,2</sub><br>(<i>M</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td>3.0</td><td>0.08</td><td>1.49</td><td>1.45</td><td>0.052</td><td>2.3e-3</td><td>1.31</td></tr><tr><td>3.0</td><td>0.10</td><td>1.58</td><td>1.50</td><td>0.063</td><td>2.0e-3</td><td>1.35</td></tr><tr><td>3.0</td><td>0.50</td><td>1.73</td><td>1.56</td><td>0.311</td><td>7.2e-4</td><td>1.40</td></tr><tr><td>3.0</td><td>2.00</td><td>1.80</td><td>1.57</td><td>1.23</td><td>4.8e-4</td><td>1.40</td></tr><tr><td>3.0</td><td>5.00</td><td>1.86</td><td>1.57</td><td>4.61</td><td>4.1e-4</td><td>1.40</td></tr><tr><td>3.0</td><td>10.0</td><td>2.00</td><td>1.57</td><td>9.51</td><td>2.3e-4</td><td>1.40</td></tr><tr><td>3.0</td><td>20.0</td><td>2.17</td><td>1.58</td><td>20.0</td><td>1.4e-4</td><td>1.41</td></tr><tr><td>3.0</td><td>50.0</td><td>2.32</td><td>1.58</td><td>53.2</td><td>6.3e-5</td><td>1.41</td></tr><tr><td>3.0</td><td>80.0</td><td>2.37</td><td>1.57</td><td>87.1</td><td>5.2e-5</td><td>1.40</td></tr><tr><td>3.0</td><td>100.0</td><td>2.38</td><td>1.57</td><td>109.7</td><td>4.8e-5</td><td>1.40</td></tr></tbody></table>

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

maximum NS accretion rate in units of the classical Eddington limit, $X_{\text{Edd}} \equiv \dot{M}_{\text{NS}}^{\text{max}}/\dot{M}_{\text{Edd}}$. The larger the value of $X_{\text{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 $\Delta M_{\text{NS}}$ in Tables 5 and 6 are based on an assumed accretion efficiency value of $X_{\text{Edd}} = 1$. However, as discussed in Section 4.4, it is possible that the amount of material accreted by the NS, $\Delta M_{\text{NS}}$, is larger by a factor of three, based on evidence for high accretion efficiency in PSR J1952+2630 (Lazarus et al. 2014).

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 ~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}_{\text{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_{\text{NS}} = M_{\text{core,f}} - M_{\text{NS}}^{\text{bind}}$, where the released binding energy of the NS corresponds to a mass decrease of $M_{\text{NS}}^{\text{bind}} = 0.084 M_{\odot} (M_{\text{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_{\text{NS,2}}$ in Table 5), given the relatively small initial helium star masses, $M_{\text{He,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_{\text{orb}}$.

At first look from Figure 8, one could be tempted to conclude that theoretical modeling with $X_{\text{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_{\text{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.

<figure><img src="image_11.png" /></figure>

**Figure 9.** Spread in the ($P_{\text{orb}}, P$) correlation (blue stars) when applying a wide range of initial binary parameters ($M_{\text{He,i}}$ and $P_{\text{orb,i}}$), which lead to quite different masses and orbital periods of the exploding helium stars. Calculations are based on data in Table 6 and assuming a purely symmetric SN ($w = 0$). The observational data are plotted with red points, and the green line is our birth relation based on this data (see Figure 6).

<table><caption>**Table 6**
Properties of Various Case BB RLO Systems</caption><thead><tr><th><i>M</i><sub>He,i</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>orb,i</sub><br>(days)</th><th><i>M</i><sub>He,f</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>core,f</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>P</i><sub>orb,f</sub><br>(days)</th><th>Δ<i>M</i><sub>NS</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>M</i><sub>NS,2</sub><br>(<i>M</i><sub>☉</sub>)</th></tr></thead><tbody><tr><td>2.8</td><td>0.10</td><td>1.43</td><td>1.39</td><td>0.079</td><td>2.7e-3</td><td>1.26</td></tr><tr><td>3.0</td><td>0.10</td><td>1.58</td><td>1.50</td><td>0.063</td><td>2.0e-3</td><td>1.35</td></tr><tr><td>3.5</td><td>0.10</td><td>1.88</td><td>1.76</td><td>0.048</td><td>1.2e-3</td><td>1.56</td></tr><tr><td>4.0</td><td>0.10</td><td>2.12</td><td>1.96</td><td>0.048</td><td>9.4e-4</td><td>1.71</td></tr><tr><td>6.0</td><td>0.10</td><td>2.37</td><td>2.15</td><td>0.064</td><td>6.7e-4</td><td>1.86</td></tr><tr><td>2.7</td><td>0.50</td><td>1.48</td><td>1.41</td><td>0.415</td><td>8.4e-4</td><td>1.27</td></tr><tr><td>3.0</td><td>0.50</td><td>1.73</td><td>1.56</td><td>0.311</td><td>7.2e-4</td><td>1.40</td></tr><tr><td>3.5</td><td>0.50</td><td>2.07</td><td>1.80</td><td>0.365</td><td>5.0e-4</td><td>1.59</td></tr><tr><td>2.6</td><td>2.00</td><td>1.46</td><td>1.37</td><td>1.78</td><td>8.3e-4</td><td>1.24</td></tr><tr><td>3.0</td><td>2.00</td><td>1.80</td><td>1.57</td><td>1.23</td><td>4.8e-4</td><td>1.40</td></tr><tr><td>3.5</td><td>2.00</td><td>2.39</td><td>1.81</td><td>1.78</td><td>2.6e-4</td><td>1.60</td></tr><tr><td>4.0</td><td>2.00</td><td>2.95</td><td>2.02</td><td>2.45</td><td>1.4e-4</td><td>1.76</td></tr><tr><td>2.6</td><td>20.0</td><td>1.56</td><td>1.37</td><td>19.3</td><td>6.3e-4</td><td>1.24</td></tr><tr><td>3.0</td><td>20.0</td><td>2.17</td><td>1.58</td><td>20.0</td><td>1.4e-4</td><td>1.41</td></tr><tr><td>3.2</td><td>20.0</td><td>2.67</td><td>1.70</td><td>23.8</td><td>2.8e-5</td><td>1.51</td></tr><tr><td>3.0</td><td>50.0</td><td>2.32</td><td>1.58</td><td>53.2</td><td>6.3e-5</td><td>1.41</td></tr><tr><td>3.0</td><td>80.0</td><td>2.37</td><td>1.57</td><td>87.1</td><td>5.2e-5</td><td>1.40</td></tr></tbody></table>

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

To illustrate the spread in theoretically estimated points in the ($P_{\text{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_{\text{He,i}} = 2.6–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–2.37 $M_{\odot}$, with metal cores of 1.37–2.15 $M_{\odot}$. The resulting NS masses cover a range between 1.2–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_{\text{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
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Eccentricity as a function of orbital period for all confirmed DNS systems observed in the Galactic disk (red points). Black dots are simulated systems based on models from Tauris et al. (2015), using $M_{\text{He,i}} = 3.0 M_{\odot}$ and $P_{\text{orb,i}} = 0.08–100$ days. For each model, 1000 resulting bound DNS systems were simulated based on SNe with a fixed kick of $w = 50$ km s$^{-1}$ and applied in a random (isotropic) direction. Unbound systems are not included. The dashed blue line is connecting discrete data point solutions for the resulting correlation in the case of symmetric SNe ($w = 0$).</figcaption></figure>

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_{\text{orb}}$ or $P$.

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

### 5.3. The ($P_{\text{orb}}$, $e$) Correlation

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_{\text{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_{\text{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_{\text{He,i}}$. Furthermore, it is more sensitive to even small kicks compared to the ($P_{\text{orb}}$, $P$) correlation.

Indeed, the observed spread in eccentricities for a given value of $P_{\text{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$ km s$^{-1}$. In Figure 10, we simulated the dynamical effects of applying kicks of $w = 50$ km 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

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Dependence of the post-SN eccentricity on the direction of the kick (see Figure 13). About half (56%) of the systems remain bound, with an average 3D systemic velocity of $\langle v_{\text{sys}} \rangle = 40$ km s$^{-1}$. The value of their eccentricities ($0 < e < 1$) are color-coded; see the scale. All systems with $\theta < \theta_{\text{crit}}$, however, will not survive the kick. The small black areas in the center of the red regions mark the almost-circular post-SN systems with $e < 0.015$. The three arrows on the left side mark (bottom to top) $\theta_{\text{crit}}$, $\theta_a$, and $\theta_p$—see Section 6. For another version of this plot, see Figure 40.</figcaption></figure>

area in the ($P_{\text{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_{\text{orb}} = 0.323$ days and $e = 0.617$) and PSR J1811−1736 ($P_{\text{orb}} = 18.8$ 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_{\text{He,i}}$, whereas PSR B1913+16 must have had $w \gtrsim 200$ km 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 < e < 1$, or even disrupt the system. This is demonstrated in Figure 11, where we modeled the explosions of $M_{\text{He,f}} = 3.0 M_{\odot}$ stars with a fixed kick magnitude of $w = 50$ km s$^{-1}$ and varied the direction of the applied kicks by systematically changing the two kick angles, $\theta$ and $\phi$ (see Section 6 for a definition). The pre-SN orbital period was assumed to be $P_{\text{orb,f}} = 25$ days in all cases. It is seen that systems with $\theta$ less than some critical value, $\theta_{\text{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).

### 5.4. The ($P$, $e$) Correlation

In the previous sections, we have investigated the correlations for DNS systems between $P_{\text{orb}}$ and $P$, and between $P_{\text{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
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Eccentricity as a function of spin period for all first-born (recycled) NSs in DNS systems. The observational data are shown with red dots. The error bars are too small to be seen. Black dots along discrete vertical lines are simulated systems based on models from Tauris et al. (2015); see Figure 10 for details. The dashed and solid blue lines show the resulting correlation for $w = 0$, and the average values of the eccentricity for bound systems using $w = 50$ km s$^{-1}$, respectively.</figcaption></figure>

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$ km s$^{-1}$. The solid blue line represents the outcome of our theoretical models with applied kicks of $w = 50$ km s$^{-1}$ averaged over all kick directions that resulted in bound systems. We assumed again an Eddington accretion efficiency parameter $X_{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_{orb,i}$, assuming a random (isotropic) kick direction for a fixed kick magnitude of $w = 50$ km 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 \le 50$ km s$^{-1}$. However, there are further observational constraints other than $P$, $P_{orb}$, and eccentricity, such as proper motion, which in some cases require larger values of $w$ (see discussions in Sections 6 and 8).

In the previous four subsections (Sections 5.1–5.4), we have discussed the correlations between $P$, $P_{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

The importance of the kick magnitude during the second SN in DNS systems for the suggested correlation between $P$ and eccentricity was investigated in detail by Dewi et al. (2005) using a population synthesis method, also taking into account the orbital and spin evolution since the birth of a given DNS system. They found that large kicks (using a Maxwellian kick velocity distribution with a 1D dispersion of $\sigma = 190$ km s$^{-1}$) would completely randomize any correlation, whereas small kicks with typically $w < 50$ km s$^{-1}$ ($\sigma = 20$ km s$^{-1}$) preserve a correlation.

In addition, Dewi et al. (2005) found that the ($P$, $e$) correlation cannot be explained simply by selection effects. Such a hypothesis had been suggested by Chaurasia & Bailes (2005), who argued that the correlation is, at least partly, caused by close-orbit (thus small $P$) systems with high eccentricities merging on a short timescale due to GWs, thereby removing them from the observed sample. This hypothesis would then explain an observed DNS population with some correlation between $P$ and eccentricity. However, using population synthesis, Dewi et al. (2005) showed that although close-orbit DNS systems are removed from the expected observed sample due to merger events, the remaining systems still show a random distribution of eccentricities for a given value of $P$ if the applied kicks are large ($\sigma = 190$ km s$^{-1}$, corresponding to average kicks of $\langle w \rangle = 450$ km s$^{-1}$). Only by applying a small kick distribution ($\sigma = 20$ km s$^{-1}$, $\langle w \rangle = 47$ km s$^{-1}$) does a correlation between $P$ and eccentricity remain, in agreement with our findings. As we shall discuss further below, a small kick imparted onto the newborn NS in the second SN is indeed, in general, providing a better fit to observations (see also Bray & Eldridge 2016). Further discussions on the selection effects of the observed population due to kicks are discussed in Section 6.9.

## 6. General Dynamical Effects of SNE

The dynamical effects of SNe in close binaries (and hierarchical triples) have been studied both analytically and numerically in a number of papers over the last four decades, since the discovery of the Hulse–Taylor pulsar, e.g., Flannery & van den Heuvel (1975), Sutantyo (1978), Hills (1983), Brandt & Podsiadlowski (1995), Tauris & Bailes (1996), Kalogera (1996), Tauris & Takens (1998), Wex et al. (2000), and Pijloo et al. (2012).

To solve for the dynamical effects of the SN explosion, the SN event can be assumed to be instantaneous given that the SN ejecta velocity is much greater than the binary orbital velocity.<fnref n="18" /> The mass loss reduces the absolute value of the potential energy and affects the orbital kinetic energy by decreasing the reduced mass of the system. In addition, a kick is imparted on the newborn NS. Considering the change in total energy of the system, the change in the orbital semimajor axis (ratio of final to initial value) can be expressed as (Hills 1983)

$$ \frac{a_f}{a_i} = \left[ \frac{1 - \Delta M / M}{1 - 2\Delta M / M - (w/v_{\text{rel}})^2 - 2\cos\theta(w/v_{\text{rel}})} \right], \quad (8) $$

where $\Delta M$ is the amount of instantaneous mass loss from the exploding star (in our notation for the second SNe applied here, $\Delta M = M_{\text{He,f}} - M_{\text{NS,2}}$), $M = M_{\text{He,f}} + M_{\text{NS,1}}$ is the total mass of the pre-SN system, and $v_{\text{rel}}$ is the relative velocity between

<footnote n="18">¹⁸ We ignore the extremely rare cases discussed in Section 4.5.</footnote>
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Illustration of the geometry of the two kick angles ($\theta$, $\phi$) and the kick sphere surrounding the exploding star (blue). The pre-SN orbital angular momentum vector is along the z-axis. The kick velocity vector ($w$) is shown in red.</figcaption></figure>

the two stars ($\sqrt{G(M_{\text{He,f}} + M_{\text{NS,1}})/a_i}$). The kick angle $\theta$ is defined as the angle between the kick velocity vector, $w$, and the pre-SN orbital velocity vector of the exploding star, $v_{\text{He}}$, in the pre-SN center-of-mass rest frame (see Figure 13). Using Kepler's third law, the change in $P_{\text{orb}}$ can be obtained.

Equation (8) applies to a circular pre-SN binary. This is most likely a good approximation here given the tidal interactions during Case BB RLO. For the more general case, see Hills (1983). While the shell-impact effects on the companion star from exploding helium stars in tight binaries can potentially be important for the first SN explosion (e.g., Wheeler et al. 1975; Liu et al. 2015; Tauris 2015), we disregard such effects in the second SN where the companion star is a NS (essentially a point mass).

Solving for the denominator being equal to zero in Equation (8) yields the critical angle, $\theta_{\text{crit}}$, so that $\theta < \theta_{\text{crit}}$ will result in the disruption of the orbit; see Figure 11. Thus, the probability of a binary system surviving a SN with a kick in a random (isotropic) orientation can be found by integration and yields $P_{\text{bound}} = 1 - (1 - \cos \theta_{\text{crit}})/2$, or equivalently (Sutantyo 1978; Hills 1983)

$$P_{\text{bound}} = \frac{1}{2} \left\{ 1 + \left[ \frac{1 - 2\Delta M/M - (w/v_{\text{rel}})^2}{2(w/v_{\text{rel}})} \right] \right\}, \quad (9)$$

where $P_{\text{bound}}$ is restricted to [0, 1] and takes the value of 0 or 1, below and above this interval, respectively.

The eccentricity of the post-SN system can be evaluated directly from the post-SN orbital angular momentum, $L_{\text{orb,f}}$, and is given by

$$e = \sqrt{1 + \frac{2 E_{\text{orb,f}} L_{\text{orb,f}}^2}{\mu_f G^2 M_{\text{NS,1}}^2 M_{\text{NS,2}}^2}}, \quad (10)$$

where

$$L_{\text{orb,f}} = a_i \mu_f \sqrt{(v_{\text{rel}} + w \cos \theta)^2 + (w \sin \theta \sin \phi)^2}. \quad (11)$$

Here, $\mu_f$ and $E_{\text{orb,f}} = -GM_{\text{NS,1}}M_{\text{NS,2}}/2a_f$ are the post-SN reduced mass and orbital energy, respectively. The kick angle $\phi$ is measured in the plane perpendicular to the pre-SN velocity vector of the exploding star, $v_{\text{He}}$ (Hills 1983; Tauris & Takens 1998), such that the component of the kick velocity pointing directly toward the companion star is given by $w_y = w \sin \theta \cos \phi$ (Figure 13).

Solving for the right-hand side in Equation (8) being equal to 1 yields another interesting angle, $\theta_a$,

$$\theta_a = \cos^{-1} \left( - \frac{\Delta M/M + (w/v_{\text{rel}})^2}{2(w/v_{\text{rel}})} \right), \quad (12)$$

which determines whether the post-SN semimajor axis, $a_f$, is larger or smaller than the pre-SN orbital radius, $a_i$. Thus, as a consequence of the SN explosion, the orbit widens<fnref n="19" /> if the kick angle $\theta < \theta_a$ and it shrinks if $\theta > \theta_a$. The ratio ($a_f/a_i$) is a simple monotonically decreasing function of $\theta$ (see Equation (8)) such that minimum-size post-SN orbits are achieved for $\theta \to 180^{\circ}$. For random kick directions, we find that the probability that the post-SN orbit shrinks ($a_f < a_i$) is given by<fnref n="20" />

$$P_{\text{a}}^{-} = \frac{1}{2} - \left[ \frac{\Delta M/M + (w/v_{\text{rel}})^2}{4(w/v_{\text{rel}})} \right]. \quad (13)$$

This expression is convenient to use as an indicator for estimating the distribution of merger timescales of DNS systems produced by the second SN.

Combining Equations (10) and (12), we can thus derive the two fine-tuned solutions for the kick angles for which the post-SN system is circular (i.e., $e=0$). The result is simply $(\theta, \phi) = (\theta_a, \pm 90^{\circ})$; see Figure 11 for an example.

As a final twist in the tail, using Kepler's third law we can also work out a critical angle, $\theta_P$,

$$\theta_P = \cos^{-1} \left( \frac{1 - 2\Delta M/M - (1 - \Delta M/M)^{2/3} - (w/v_{\text{rel}})^2}{2(w/v_{\text{rel}})} \right), \quad (14)$$

above which the post-SN orbital period is smaller than the pre-SN orbital period. The probability that $\theta > \theta_P$ (and therefore that $P_{\text{orb}}$ decreases) is thus given by

$$P_{\text{P}}^{-} = \frac{1}{2} - \left[ \frac{2\Delta M/M + (1 - \Delta M/M)^{2/3} + (w/v_{\text{rel}})^2 - 1}{4(w/v_{\text{rel}})} \right]. \quad (15)$$

It is always the case that $\theta_{\text{crit}} < \theta_a < \theta_P$ for any value of $w$ and $\Delta M > 0$. Thus, for systems with a kick angle $\theta_a < \theta < \theta_P$, the orbital semimajor axis shrinks, while at the same time, the orbital period increases as a result of the SN. Examples of $\theta_{\text{crit}}$, $\theta_a$, and $\theta_P$ for a given SN explosion are shown in Figure 11. Finally,

<footnote n="19">¹⁹ It should be noted that in case the term in the parentheses in Equation (12) takes a value less than -1, it follows that $\theta_a = 180^{\circ}$, meaning that all post-SN systems will widen.</footnote>

<footnote n="20">²⁰ $P_a^-$ is restricted to $[0, \frac{1}{2}]$ and is always 0 below this interval.</footnote>
constraints on retrograde versus prograde post-SN spin orbits are discussed in Hills (1983) and Brandt & Podsiadlowski (1995).

## 6.1. Misalignment Angles

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 \neq 0^{\circ}$ and $\phi \neq \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)

$$ \delta = \tan^{-1} \left( \frac{w \sin \theta \sin \phi}{\sqrt{(v_{\text{rel}} + w \cos \theta)^2 + (w \sin \theta \cos \phi)^2}} \right) \quad (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 < 3^{\circ}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 $\gamma$-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 $\delta$ of the recycled NSs to constrain kicks in the second SN event.

## 6.2. Systemic Velocities

Another important diagnostic quantity (aside from $P_{\text{orb}}$, P, eccentricity, $M_{\text{NS}}$, and $\delta$) for understanding DNS formation is their post-SN systemic velocity, $v_{\text{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

$$ v_{\text{sys}} = \sqrt{(\Delta p_x)^2 + (\Delta p_y)^2 + (\Delta p_z)^2} / (M_{\text{NS},1} + M_{\text{NS},2}), \quad (17) $$

where the change in momentum is given by

$$ \begin{aligned} \Delta p_x &= M_{\text{NS},2} w \cos \theta - \Delta M M_{\text{NS},1} \sqrt{G/(Ma_i)}, \\ \Delta p_y &= M_{\text{NS},2} w \sin \theta \cos \phi, \\ \Delta p_z &= M_{\text{NS},2} w \sin \theta \sin \phi. \end{aligned} \quad (18) $$

<figure><img src="image_16.png" /><figcaption>Figure 14. Distribution of $v_{\text{sys}}$ for the simulated DNS systems shown in Figure 10 (using $w = 50 \text{ km s}^{-1}$); see the text for discussions.</figcaption></figure>

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 \text{ km s}^{-1}$ were applied. For these resulting systems, the average systemic velocity is $\langle v_{\text{sys}} \rangle \approx 30 \text{ km 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 < 50 \text{ km s}^{-1}$ (see detailed discussions further below).

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-20 \text{ km 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).

## 6.3. Symmetric SNe

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)

$$ \frac{a_f}{a_i} = \left[ \frac{1 - \Delta M / M}{1 - 2 \Delta M / M} \right] \quad (19) $$

and

$$ e = \Delta M / (M_{\text{NS},1} + M_{\text{NS},2}), \quad (20) $$

and where the probability of remaining bound is always $P_{\text{bound}} = 1$ for $\Delta M / M < 0.5$, whereas (following the virial theorem) all systems are disrupted if more than half of the total mass is lost, i.e., if $\Delta M / M > 0.5$.
### 6.3.1. The Eccentricity Floor and the NS EoS

If the second SN in DNS systems were always symmetric ($w=0$) and the amount of SN ejecta mass negligible (i.e., assuming extremely ultra-stripped exploding cores), then the only contribution to $\Delta M$ would be the released gravitational binding energy ($M_{NS}^{bind}$). This is released when baryonic mass is converted into gravitational mass during the collapse (see, e.g., Section 5.2). Under such circumstances, one would expect to observe an “eccentricity floor” at approximately $e = \Delta M/M \simeq 0.06$. From observations of many DNS systems, the exact inferred value of this “eccentricity floor” would provide a measure of the released gravitational binding energy, $E_{bind}$, and thereby constrain the NS EoS (Lazarus et al. 2016).

Unfortunately, this method is difficult in practice since even ultra-stripped NS progenitors have small amounts ($<0.20 M_{\odot}$) of baryonic SN ejecta mass. For symmetric explosions with such small amounts of ejecta mass, the post-SN eccentricities of DNS systems are expected to be $e \simeq 0.07-0.13$, in agreement with constraints obtained from some of the known DNS systems (see Section 8). Such low eccentricities are, however, also possible for asymmetric SNe, depending on the kick direction. Furthermore, even small kicks of $w \sim 50$ km s$^{-1}$ will generally lead to a substantial spread in eccentricities (also to cases with $e < 0.06$), as we have demonstrated here in Figure 10, thereby making the idea of constraining the NS EoS impossible in practice using this method. Finally, in sufficiently tight binaries, gravitational damping may have decreased the values of the eccentricity significantly since the second SN explosion (Lazarus et al. 2016).

### 6.4. NS Kick Magnitudes in Binaries

Whereas average kicks of 400–500 km 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–500 km 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$ km 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.

### 6.5. Theoretical Estimates of NS Kick Magnitudes

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

$$ \xi_M = \frac{M/M_{\odot}}{R(M)/1000 \text{ km}}, \qquad (21) $$

where $R(M)$ is the radius of the coordinate enclosing mass, $M$. Sukhbold & Woosley (2014) argued that $\xi$ at bounce is correlated with $\xi$ 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, $\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_{\text{exp}}$:

$$w = 211 \text{ km s}^{-1} C \left(\frac{\alpha_{\text{ej}}}{0.1}\right) \left(\frac{E_{\text{exp}}}{10^{51} \text{ erg}}\right) \left(\frac{M_{\text{NS}}}{1.5 M_{\odot}}\right)^{-1}, \quad (22)$$

where C is a constant typically of order unity, $M_{\text{NS}}$ is the (baryonic) NS mass, and $\alpha_{\text{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_{\text{exp}} \simeq 3 \times 10^{50}-2.5 \times 10^{51} \text{ erg}$ and $\alpha_{\text{ej}} = 0-0.33$, yielding kick values in a broad range between $w = 0-1000 \text{ km 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_{\text{exp}} \simeq 10^{50} \text{ erg}$ and $\alpha_{\text{ej}} \lesssim 0.03$, resulting in very small kicks of $w \lesssim 10 \text{ km s}^{-1}$.

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

The concept of ultra-stripped SNe was recently addressed in Tauris et al. (2013, 2015). These are SNe where a naked helium star experiences extreme stripping by a close-orbit compact object prior to its explosion. Ultra-stripped SNe are therefore particularly relevant for the second SN in DNS systems (Figure 1). Their expected optical signatures are compatible with observations of rapidly decaying SN light curves, such as SN 2005ek (Drout et al. 2013; Tauris et al. 2013; Moriya et al. 2017).

The flavor of ultra-stripped SNe can be either an EC SN, from a collapsing ONeMg core (Nomoto 1987), or an Fe CCSN. The momentum kick imparted on a newborn NS via an EC SN is always expected to be small. This follows from detailed simulations, which imply (i) explosion energies significantly

<figure><img src="image_17.png" /><figcaption>Figure 15. Dependence of the estimated NS kick magnitude, w, on the SN explosion energy, $E_{\text{exp}}$, and the momentum-asymmetry parameter, $\alpha_{\text{ej}}$; see Equation (22). Hyperbolas for constant values of w = 50, 200, 500, 1000, and 1500 km s$^{-1}$ are shown.</figcaption></figure>

smaller than those inferred for standard Fe CCSNe (Dessart et al. 2006; Kitaura et al. 2006), and (ii) short timescales to revive the stalled SN shock, compared to the timescales of the non-radial hydrodynamic instabilities, which are required to produce strong anisotropies; see, e.g., Podsiadlowski et al. (2004), and Janka (2012). Simulations of EC SNe by Kitaura et al. (2006) and Dessart et al. (2006) predict explosion energies of about (or even less than) $10^{50} \text{ erg}$, and thus most likely kick velocities below $50 \text{ km s}^{-1}$ (Figure 15).

Whereas Fe CCSNe are certainly able to produce large kicks (e.g., Janka 2017), small NS kicks have been suggested to originate from CCSNe with small iron cores (Podsiadlowski et al. 2004 and references therein). For such small Fe CCSNe, the situation is similar to that of EC SNe. In both cases, SN simulations (A. Gessner & H.-T. Janka 2017, in preparation) suggest fast explosions where non-radial hydrodynamical instabilities (convectively driven or from the standing accretion shock) are unable to grow, leading to somewhat small kick velocities.

Using binary stellar evolution arguments, Tauris et al. (2015) identified two factors that also imply that in ultra-stripped SNe, the NS kicks may by small. First, from their modeling of the progenitor stars of ultra-stripped SNe, they have demonstrated that the expected amount of ejecta is extremely small ($\sim 0.1 M_{\odot}$) compared to standard SN explosions in which several solar masses of material is ejected. This may likely lead to a weaker gravitational tug on the proto-NS (caused by the asphericity of the ejecta; see, e.g. Janka 2012, 2017) and thus a small kick. Second, the binding energies of the envelopes of their final progenitor star models are often only a few $10^{49} \text{ erg}$, such that even a weak outgoing shock can quickly lead to their ejection, potentially before large anisotropies can build up.

To mimic the outcome, but avoiding detailed binary stellar evolution calculations of Case BB RLO, Suwa et al. (2015) performed axisymmetric hydrodynamical simulations of neutrino-driven explosions of ultra-stripped Fe CCSNe using the stellar evolution outcomes of single, evolved CO stars. All their models exhibited successful explosions driven by neutrino heating. Their diagnostic explosion energy, ejecta mass, and nickel mass were typically $10^{50} \text{ erg}$, $0.1 M_{\odot}$, and $0.01 M_{\odot}$, respectively, i.e.,
compatible with observations of rapidly decaying light curves such as SN 2005ek (Drout et al. 2013; Tauris et al. 2013). Moreover, their calculated kick velocities were typically less than 50 km s⁻¹, and sometimes even below ~10 km s⁻¹, in agreement with the simulations discussed in Janka (2017) and references therein.

## 6.6. Relation Between NS Mass and SN Kick?

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–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

While the relation in Equation (22) is supported by numerical SN simulations (see references discussed by Janka 2017), a dependence of the NS kick on the total SN ejecta mass or NS mass is more indirect and enters through possible correlations of these quantities with $E_{\text{exp}}$. Janka (2017) provides arguments on why collapse events of stellar cores with low compactness (which are associated with EC SNe, low-mass Fe CCSNe, and ultra-stripped SNe with small metal cores) can be expected to produce NSs with considerably smaller kicks than SNe from more ordinary iron cores. Since stellar cores with low compactness also typically produce lower-mass NSs (see Ugliano et al. 2012; Ertl et al. 2016; Sukhbold et al. 2016), a relation between SN kick and NS mass seems plausible. We note that although in Equation (22) $w$ appears to scale inversely proportionally to $M_{\text{NS}}$, the values of $E_{\text{exp}}$ and $\alpha_{\text{ej}}$ are systematically smaller for explosions leading to small values of $M_{\text{NS}}$, such that a correlation between $w$ and $M_{\text{NS}}$ is indeed expected in general.

### 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–1400$ km s⁻¹, although with a peak at a low velocity of $w \simeq 70$ km s⁻¹). For the second-born NSs with smaller masses, we note that the inferred kicks are quite small (or could potentially have been small).

Our hypothesis of a possible correlation between $M_{\text{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

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Kick velocity as a function of NS mass for the second SN explosion in the observed DNS systems (Table 7). For each DNS system, we applied two distance models; see Section 2.2. For clarity, the pink triangles are plotted systematically with an artificial mass offset of +0.002 $M_\odot$. The more massive NSs (PSR B1534+12 and B1913+16) seem to have received significantly larger kicks compared to the low-mass NS systems (see Sections 6.6, 8.4, and 8.12).</figcaption></figure>

<table><caption>**Table 7**
Young NS Masses and Resulting NS Kicks in DNS Systems</caption><thead><tr><th>Pulsar</th><th><i>M</i><sub>NS,2</sub><br>(<i>M</i><sub>☉</sub>)</th><th><i>w</i><sub>peak</sub><br>(km s<sup>-1</sup>)</th><th><i>w</i><sub>interval</sub><br>(km s<sup>-1</sup>)</th><th><i>e</i></th></tr></thead><tbody><tr><td>B1913+16</td><td>1.389</td><td>300</td><td>185–465</td><td>0.617</td></tr><tr><td>B1534+12</td><td>1.346</td><td>240</td><td>175–300</td><td>0.274</td></tr><tr><td>J1906+0746</td><td>1.291</td><td>70</td><td>0–1500</td><td>0.085</td></tr><tr><td>J0737–3039</td><td>1.249</td><td>18</td><td>0–70</td><td>0.088</td></tr><tr><td>J1756–2251</td><td>1.230</td><td>10</td><td>0–130</td><td>0.181</td></tr><tr><td>J0453+1559</td><td>1.174</td><td>45</td><td>15–175</td><td>0.113</td></tr></tbody></table>

**Note.** DNS systems listed in order of measured masses of the second-born (young) NSs. The values of $w_{\text{peak}}$ (third column) denote the peak values of the solutions to their inferred kick velocity distribution with an interval of $w_{\text{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.

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
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Eccentricity as a function of second-born NS mass; see Table 7. PSR J1518+4904 is marked in magenta color, due to an uncertainty in its mass (G. Janssen 2017, private communication). Although a correlation may seem to be present in the current data, we expect a large scattering in this diagram as more DNS sources are discovered.</figcaption></figure>

state of the stellar progenitor immediately prior to the second SN event.

## 6.7. Relation Between NS Mass and Eccentricity?

Given our hypothesis of a mass-kick correlation for the second SN in DNS systems, one might also expect an ($M_{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_{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.

## 6.8. Nature of Ultra-stripped SNe in DNS Progenitors

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–11 $M_\odot$. The corresponding mass window for producing EC SNe from helium stars in close binaries is restricted to a width of ~0.2 $M_\odot$ (Tauris et al. 2015) and only for helium stars with a mass of $2.6 \le M_{He,i}/M_\odot \le 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 ($\sim 1.37–1.43 M_\odot$) do actually produce Fe CCSNe. The initial helium star mass range for producing Fe CCSNe possibly extends up to ~8–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.

## 6.9. Selection Effects from Kicks in the DNS Population

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_{orb} = 0.1$ days, the probability of surviving an ultra-stripped SN remains at 100% all the way up to $w \sim 220 \text{ km s}^{-1}$, whereas for a wide pre-SN system with $P_{orb} = 10$ days, the survival probability already declines from 100% at $w \sim 50 \text{ km s}^{-1}$, and reaches zero, with no systems surviving, for $w > 330 \text{ km s}^{-1}$. The survival probability for a more massive exploding star ($M_{He,f} = 3.0 M_\odot$, dotted line) is generally smaller because more mass is lost in the SN—this amount ($\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.

Since the majority of the observed DNS systems seem to have had a pre-SN orbital period of 0.1–1 days (based on the observed distribution of $P_{orb}$ for DNS systems and our simulations in Section 8), the survival probability remains within roughly 50% even for large kicks up to $400 \text{ km 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 \text{ km 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_{gwr}$ (Section 7), it is often useful to know the probability of the orbit shrinking as a consequence of the SN, i.e., $(a_f/a_i) < 1$, leading to small values of $\tau_{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_{orb} = 0.1$ days, it is obvious that for small kick
<figure><img src="image_20.png" /><figcaption>Figure 18. Probability of surviving the second SN in a binary as a function of the imparted NS kick velocity, w, and for three different pre-SN orbital periods of 0.1, 1, and 10 days—see Equation (9). The solid and dashed lines are for ultra-stripped SNe with ($M_{He,f}/M_{\odot}$, $\Delta M/M_{\odot}$) = (1.50, 0.30) and (1.80, 0.45), respectively, and the dotted line is for a less stripped SN (3.00, 1.45). In all cases, we assumed a first-born NS mass of $M_{NS,1} = 1.40 M_{\odot}$.</figcaption></figure>

values ($w \lesssim 50$ km s⁻¹) the post-SN systems will always become wider ($P_a = 0$). For very wide-orbit pre-SN systems, e.g., $P_{orb} = 10$ days, it is the large kicks (here $w \simeq 280-330$ km s⁻¹) 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 km s⁻¹, the pre-SN systems with $P_{orb} = 10$ days are all disrupted as shown in Figure 18.

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_a^-$ can easily be derived from Equation (13). We find that the location of the peak is at ($w/v_{rel}$) = $\sqrt{\Delta M/M}$, and its value is given by

$$P_a^-(\text{max}) = \frac{1}{2} \left[ 1 - \sqrt{\frac{\Delta M}{M}} \right]. \qquad (23)$$

As an example, for the solid lines in Figure 19, we have ($M_{He,f}/M_{\odot}$, $\Delta M/M_{\odot}$) = (1.50, 0.30) and $M_{NS,1} = 1.40 M_{\odot}$, which yields $P_a^-(\text{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 $\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 $\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_{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

<figure><img src="image_21.png" /><figcaption>Figure 19. Probability of the post-SN orbit shrinking as a consequence of the SN explosion with a kick value w (i.e., ($a_f/a_i$) < 1; see Equations (8) and (13)) for the same systems plotted in Figure 18.</figcaption></figure>

remain bound (and later evolve to HMXB sources). The solid, dashed, and dotted lines are for different SNe Ib with typical values of ($M_{He,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_{orb} = 120$ days will survive a kick of $w > 300$ km s⁻¹. This selection effect must be kept in mind when trying to reconstruct the kick magnitudes based on observations of HMXB systems.

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).

## 7. Mapping Observed DNS Systems to Simulated Post-SN DNS Systems

When comparing observational data of DNS systems with theoretically derived properties, it is important to take into account selection effects. As discussed already in Section 2.1, radio pulsars in highly relativistic (tight) orbit DNS systems with $P_{orb}$ less than a few hours are more difficult to detect than radio pulsars in wider orbits. Furthermore, as pointed out
<figure><img src="image_22.png" /><figcaption>**Figure 20.** Probability of surviving the first SN in a binary as a function of the imparted NS kick velocity, $w$, and for three different pre-SN orbital periods of 3, 30, and 120 days—see Equation (9). The solid, dashed, and dotted lines are for SNe Ib with ($M_{He,f}/M_{\odot}$, $M_2/M_{\odot}$) = (4.0, 22.0), (2.0, 15.0), and (4.0, 12.0), respectively. We assumed a newborn NS mass of $M_{NS,1} = 1.40 M_{\odot}$.</figcaption></figure>

above, the DNS systems in close orbits evolve with time. This means that the post-SN parameters applicable to DNS systems at birth ($P$, $P_{orb}$, and eccentricity) are often somewhat different from their current observed values.

As an example, consider PSR B1913+16 (the Hulse–Taylor pulsar). In Figure 21, we have plotted the past and the future evolution of this system with respect to its semimajor axis and orbital eccentricity using the well-known quadrupole formalism of GR (Einstein 1918; Peters 1964; Landau & Lifshitz 1971; Blanchet 2014). This lowest-order secular orbital evolution is easily calculated from changes in the elements of the relative orbit of two point masses resulting from GW damping, and the rates of change are (Peters 1964)

$$ \left\langle \frac{da}{dt} \right\rangle = -\frac{64}{5} \frac{G^3 M^2 \mu}{c^5 a^3 (1 - e^2)^{7/2}} \left( 1 + \frac{73}{24} e^2 + \frac{37}{96} e^4 \right), \quad (24) $$

$$ \left\langle \frac{de}{dt} \right\rangle = -\frac{304}{15} \frac{G^3 M^2 \mu e}{c^5 a^4 (1 - e^2)^{5/2}} \left( 1 + \frac{121}{304} e^2 \right), \quad (25) $$

which can be transformed into an expression for the merger time, $\tau_{gwr}$, as a function of the initial values:

$$ \tau_{gwr}(a_0, e_0) = \frac{12 C_0^4}{19 \beta} \times \int_0^{e_0} \frac{e^{29/19} [1 + (121/304)e^2]^{1181/2299}}{(1 - e^2)^{3/2}} de, \quad (26) $$

where the constants are given by

$$ C_0 = \frac{a_0(1 - e_0^2)}{e_0^{12/19}} [1 + (121/304)e_0^2]^{-870/2299}, \quad (27) $$

and

$$ \beta = \frac{64G^3}{5c^5} M^2 \mu. \quad (28) $$

Equation (26) cannot be solved analytically and must be evaluated numerically, unless we consider circular orbits, which can easily be solved analytically: $\tau_{gwr}^{\text{circ}} = a_0^4/4\beta$. The

<figure><img src="image_23.png" /><figcaption>**Figure 21.** Past ($t < 0$) and future ($t > 0$) evolution of the semimajor axis (top panel) and eccentricity (bottom panel) of PSR B1913+16. The current system is shown with a solid star, located at $t = 0$. Constraints on the spin evolution (Figure 22), however, yields a maximum age of this system of about 217 Myr (marked with an arrow). The system will merge in 301 Myr.</figcaption></figure>

merger timescale is very dependent on both $a_0$ and $e_0$. Tight and/or eccentric orbits spiral-in much faster than wider and more circular orbits

The current location of PSR B1913+16 is plotted in Figure 21 by a solid star at $t = 0$ ($a_0 = 2.80 R_{\odot}$, $e_0 = 0.617$). The system is doomed to merge in 301 Myr and will become a strong Galactic GW source at that time. For comparison, the estimated Galactic merger rate of DNS systems is of the order $1–100 \text{ Myr}^{-1}$ (Voss & Tauris 2003; Abadie et al. 2010). If this system had formed 3 Gyr ago, the semimajor axis and eccentricity at birth would have been about 21 $R_{\odot}$ and 0.95, respectively—i.e., considerably larger than observed today. The question is, how can we place limits on the age of this system to learn about its post-SN birth properties and constrain the SN explosion physics. To answer this question, we need to look at the spin evolution of the recycled NS.

For magnetodipole spin-down of a pulsar (or spin-down due to plasma currents in the magnetosphere and a pulsar wind), the braking torque $N \propto B^2$ (Manchester & Taylor 1977; Lorimer & Kramer 2004). In more general terms, one can express the braking torque as $N = I \dot{\Omega} \propto \Omega^n$, where $n$ is the braking index and $\Omega = 2\pi/P$. By integrating the deceleration law, $\dot{\Omega} = -K\Omega^n$, one obtains the following solution at time $t$ (positive in the future;
<figure><img src="image_24.png" /><figcaption>**Figure 22.** Past and future spin evolution of the recycled pulsar in PSR B1913+16, calculated for three different values of the braking index, $n = \{2, 3, 5\}$. The current system is shown with a solid star, located at $t = 0$. The system will merge in 301 Myr.</figcaption></figure>

negative in the past, e.g., Lazarus et al. 2014):

$$P = P_0 \left[ 1 + (n-1) \frac{\dot{P}_0}{P_0} t \right]^{1/(n-1)}, \quad (29)$$

$$\dot{P} = \dot{P}_0 \left( \frac{P}{P_0} \right)^{2-n}. \quad (30)$$

In Figure 22, we have plotted the past and the future spin evolution of PSR B1913+16 using three different values of the braking index, $n = \{2, 3, 5\}$; see, e.g., Tauris et al. (2012) and references therein for a discussion. Given the relatively large B-field of this pulsar ($B_0 \simeq 7 \times 10^9$ G), the braking torque acting on it is substantial. Hence, even in the most conservative case, assuming evolution with a constant $\dot{P}$ ($n = 2$), the maximum age of PSR B1913+16 is less than about 217 Myr (twice the characteristic age; see Equation (29)), given that the spin period after recycling must have been larger than zero. Hence, returning to Figure 21 (see the arrows), we find that the maximum values of the semimajor axis and the eccentricity after the SN in this case would be $a = 3.88 R_\odot$ (corresponding to $P_{\text{orb}} = 12.6$ hr) and $e = 0.710$, respectively. However, there is little evidence for pulsars generally evolving with a such a small long-term braking index of $n = 2$. Tauris & Konar (2001) argued in favor of long-term evolution with $n > 3$ (see further support in Johnston & Karastergiou 2017) and showed as an example that a pulsar born with $P_0 = 10$ ms and $\dot{P}_0 = 10^{-12.5} \text{ s s}^{-1}$ would spin down to $P = 10$ s in only 1 Myr, in case it evolved with a constant $n = 2$ (see also Equation (29)). This is clearly not the case for radio pulsars in general, and therefore we shall assume a standard evolution with $n = 3$. In this case, the maximum age of the system since recycling is given exactly by the characteristic age of the recycled pulsar (see Equation (29) for $n = 3$ and solve for $P = 0$). For PSR B1913+16, this yields a maximum age of $t = \tau_0 = 108.7$ Myr and thus upper limits on the post-SN parameters of $a = 3.34 R_\odot$ (corresponding to $P_{\text{orb}} = 10.1$ hr) and $e = 0.670$. In the following section, we take into account this

uncertainty in the age of a given DNS system when analyzing its SN explosion properties.

Finally, for each DNS system, we need to relate the observed projected 2D velocity ($v_\perp = v_T^{\text{SBB}}$) in the plane of the sky to its underlying 3D systemic velocity ($v_{\text{sys}} = v^{\text{LSR}}$). In general, assuming an isotropic orientation of $v_{\text{sys}}$ and given $v_\perp$, we find from a simple integration that the probability that $v_{\text{sys}}$ is less than a certain value $v'_{\text{sys}}$ is given by

$$P(v_{\text{sys}} < v'_{\text{sys}}) = 1 - \sqrt{1 - \psi^2}, \quad (31)$$

where $\psi$ is the ratio $v_\perp/v'_{\text{sys}}$. However, rather than the 3D systemic velocity of the DNS system relative to the SSB, we are interested in the 3D systemic velocity in its LSR, $v^{\text{LSR}}$. These velocities were calculated via the method explained in Section 2.2. In the simulations below, we use the derived velocities ($v^{\text{LSR}}$) from Table 3 to help constrain the SN explosion and the pre-SN binary parameters.

## 8. Simulations of the Second SN in DNS Progenitor Systems

We apply MC techniques to simulate SN explosions (between 10 and 500 million trials for each system) to calculate the kinematic properties of DNS systems surviving the second SN explosion. From the outcome of these simulations, we can compare with any given observed DNS system and iterate on adjusting the pre-SN parameter space until the outcome matches the observed post-SN values within a certain error margin (including the evolution of $P_{\text{orb}}$ and eccentricity since the formation of the DNS system, if necessary; see Section 7).

Our simulations take their basis in a five-dimensional phase space. The input parameters are the pre-SN orbital period, $P_{\text{orb},1}$; the final mass of the (stripped) exploding star, $M_{\text{He,f}}$; the magnitude of the kick velocity imparted on the newborn NS, $w$; and the two angles defining the direction of the kick velocity, $\theta$ and $\phi$. A sixth input parameter is the mass of the first-born NS, $M_{\text{NS},1}$. However, the SN simulation results are not very dependent on the small variations of this parameter.

For almost all DNS systems, the total mass, $M$, is known from measurements of the periastron advance of the radio pulsar:

$$\dot{\omega} = 3 T_{\odot}^{2/3} \left( \frac{P_{\text{orb}}}{2 \pi} \right)^{-5/3} \frac{M^{2/3}}{1 - e^2}, \quad (32)$$

assuming that this effect is caused purely by the effects of GR to first post-Newtonian order (Robertson 1938; Taylor & Weisberg 1982), and where $T_\odot \equiv GM_\odot c^{-3} = 4.925490947 \text{ µs}$. Hence, we can use this post-SN total mass to constrain the combination of $M_{\text{NS},1}$ and $M_{\text{He,f}}$ prior to the SN explosion by requiring that $M_{\text{He,f}}$ cannot be less than $M_{\text{NS},2}$ (after correction for the released gravitational binding energy; see Section 5.2).

To compute the outcome of an asymmetric SN in a binary system, analytic expressions for the general solutions have been found both for systems that remain bound (Hills 1983) and for systems that disrupt (Tauris & Takens 1998). We follow the notation and the methodology outlined in these papers as summarized in Section 6. The formulae of Tauris &
Takens (1998) are solutions for the general case and reproduce the formulae of Hills (1983) for bound systems.

In Figures 24–38 in the Appendix, we show the results of our simulations of each observed DNS system in six panels. The upper-left panel shows the post-SN values of $P_{\text{orb}}$ and eccentricity. The chosen solutions are plotted as small black dots (making up a square) centered on the observed values of ($P_{\text{orb}}$, $e$) with an accepted error margin of ±3% in both $P_{\text{orb}}$ and $e$. All constraints on the observed parameters of the DNS system are also written in this panel, as well as the number of simulated SN explosions, $N$.

The upper-right panel shows the distribution of the resulting 3D systemic velocities. These velocities can be combined with the merger time for the DNS system to estimate the offset distance of the future DNS merger (short GRB) from its birth place in the Galactic disk (e.g., Bloom et al. 1999; Perna & Belczynski 2002; Voss & Tauris 2003; Fong et al. 2010; Kelley et al. 2010; Beniamini et al. 2016b). Assuming, as a first-order approximation, constant motion throughout a galactic potential, an upper limit for this distance can be simply estimated as $d_{\text{pc}} \simeq v_{\text{km/s}} \cdot \tau_{\text{gwr,Myr}}$. The recent discovery of significant abundances of r-process elements (e.g., Eu) in an ultra-faint dwarf galaxy (Reticulum II; Roederer et al. 2016) requires very small kick velocities of $w < 50 \text{ km s}^{-1}$ for the second SN explosion (Beniamini et al. 2016a), if these elements are produced by a DNS merger retained in such a small galaxy. For ultra-faint dwarf-like galaxies, the escape velocity can be quite small (down to ~15 km s$^{-1}$). However, a number of the known Galactic DNS systems could, in principle, have formed with such small 3D systemic velocities; see our simulations discussed below and the plots in Figures 24–38.

The central panels show the pre-SN orbital period, $P_{\text{orb,i}}$ (left), and the mass of the exploding star, $M_{\text{He,f}}$ (right). The intervals of the trial values are also stated in the legend, and these intervals must be broader than (and cover the entire range of) the solution values. The mass of the exploding star, $M_{\text{He,f}}$ is typically restricted between 1.50–7.0 $M_{\odot}$. The lower limit is determined by observational constraints on $M_{\text{NS,2}}$, taking the released gravitational binding energy from the core collapse into account. We assumed a conservative upper limit of 7.0 $M_{\odot}$ for $M_{\text{He,f}}$, since naked metal cores that were stripped by winds and RLO and still have a final mass above this limit will form a BH rather than a NS. It is, however, quite possible that BHs are also formed from less massive exploding stars.

The lower-left panel shows the applied kick velocities for all successful solutions, and the lower-right panel shows the two associated kick angles. The red circles indicate kick directions parallel ($\phi = 90^\circ, \theta = 90^\circ$) and anti-parallel ($\phi = -90^\circ, \theta = 90^\circ$) to the pre-SN orbital angular momentum vector; see Section 9.2.

We now proceed by discussing the solutions for all of the relevant SN parameters for each of the known DNS systems listed in Table 2. The corresponding plots are shown in the Appendix.

### 8.1. PSR J0453+1559

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 \text{ km s}^{-1}$ during the second SN (for an assumed upper limit on the resulting systemic velocity, $v_{\text{sys}} < 85 \text{ km s}^{-1}$; Table 3), although the tail of solutions extends to relatively large kick values up to $w \simeq 175 \text{ km s}^{-1}$. For the mass of the exploding star (producing the second-born NS), we find solutions for $M_{\text{He,f}} < 3.2 M_{\odot}$. The pre-SN orbital period is found to be $2.5 \lesssim P_{\text{orb,i}}/\text{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 < 100 \text{ km 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 < 75 \text{ km s}^{-1}$. An upper limit for the mass of the exploding star is in this case $M_{\text{He,f}} < 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_{\text{He,f}}$. Any corrections for gravitational damping in this system are negligible due to its wide orbit ($P_{\text{orb}} = 4.07 \text{ days}$).

### 8.2. PSR J0737−3039

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 < 3^\circ 2$ (Ferdman et al. 2013), and from VLBI observations (Deller et al. 2009) its distance ($\sim 1150 \pm 200 \text{ pc}$) and proper motion ($4.37 \pm 0.65 \text{ mas yr}^{-1}$) have been revealed. These measurements result in a small 2D transverse velocity of $\sim 25 \text{ km 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} \text{ km s}^{-1}$ (Table 3).

The resulting kick velocity found from our simulated solutions, where we restrict ourselves to only select systems with $v_{\text{sys}} = 18–57 \text{ km s}^{-1}$, is therefore $w \lesssim 70 \text{ km 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 \text{ km 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–20 km s$^{-1}$; see Section 6.2. A close inspection of Figure 25 shows that there are no solutions for $v_{\text{sys}} \le 32 \text{ km s}^{-1}$. The reason for this is the finite amount of mass loss, $\Delta M$.

For a comparison, we plot in Figure 26 our simulated solutions without applying any constraints on $v_{\text{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_{\text{He,f}}$.

Further important conclusions on the formation of PSR J0737−3039 can immediately be drawn from our simulations. Applying the constraint $v_{\text{sys}} < 57 \text{ km s}^{-1}$ provides