© ESO 2019

1. Introduction

Stars with initial mass larger than about 7.5 M_⊙ are the progenitors of black holes (BH) and neutron stars (NS). These stars play an important role in shaping galaxies through their radiative, chemical, and mechanical feedback (e.g., Larson 1974; Ceverino & Klypin 2009). Most young, unevolved massive stars have a nearby companion with which they form a close binary system (e.g., Sana et al. 2012; Chini et al. 2012; Kobulnicky et al. 2014; Almeida et al. 2017). Binary systems that remain bound throughout the entire evolution of both stars can give rise to many exotic phenomena, including X-ray binaries (e.g., Gott 1971; Bolton 1972; Webster & Murdin 1972; van den Heuvel & Heise 1972), binary neutron stars (e.g., Hulse & Taylor 1975; Wijers et al. 1992), gamma-ray bursts (e.g., Izzard et al. 2004a; Becerra et al. 2016; Abbott et al. 2017a), and gravitational wave events (Abbott et al. 2016, 2017b). However, only a small fraction of massive stars born in binary systems are expected to stay together their entire lives. The majority of systems are disrupted by the first core collapse event, which can separate the newly formed compact object from its former companion star (e.g., Tauris & Takens 1998; Belczyński & Bulik 1999; Belczynski et al. 2008; Eldridge et al. 2011).

Single NSs can be detected as pulsars or magnetars. Many of them are observed to have large proper motions (e.g., Gunn & Ostriker 1970; Lyne & Lorimer 1994; Hobbs et al. 2005, but see also Verbunt & Cator 2017). Conversely, single BHs are only detectable under special circumstances, for example through lensing events when passing in front of a background star (e.g., Wyrzykowski et al. 2016), or if they accrete gas from the ambient medium (e.g., Fender et al. 2013; Gaggero et al. 2017). Therefore, to probe the population of stellar-mass black holes we are limited to X-ray or gravitational wave observations that can only target the rare cases that remain bound to their companion and are close enough for the BH to accrete or merge. To learn more about the black holes that form in less special cases, it is worth investigating the imprints they may leave on their former companion star.

The primary focus of this study is the population of unbound main sequence companions, which can be identified observationally because of their peculiar spatial velocities compared to the surrounding population. Their velocities can either be detected as proper motions (i.e., their motion in the plane of the sky measured directly from the displacement of the star) or as radial velocities (i.e., their motion perpendicular to the plane of the sky, measured from the Doppler shift of the spectral lines). Large spatial velocities have been inferred for a significant sub-population of young massive stars (e.g., Blaauw 1961; Cruz-González et al. 1974; Gies & Bolton 1986; Gies 1987; Hoogerwerf et al. 2000, 2001; Tetzlaff et al. 2011; Boubert & Evans 2018).

Blaauw (1961) introduced the term “runaway stars” for those in the fast tail of the velocity distribution for a given spectral type. The typical threshold adopted to define the tail of this distribution for O and B-type stars is v ≳ 30 km s⁻¹ (e.g., Blaauw 1956; Gies & Bolton 1986; De Donder et al. 1997; Hoogerwerf et al. 2000, 2001; Dray et al. 2005; Eldridge et al. 2011), although sometimes other values have been considered (e.g., 40 km s⁻¹ in Blaauw 1961; de Wit et al. 2005; Boubert & Evans 2018). However, as we will argue based on simulations presented in this work, the majority of unbound companions are expected to exhibit velocities well below these thresholds. We will refer to these slow unbound former companions resulting from disrupted binary systems as “walkaway stars” (to our knowledge first coined by de Mink et al. 2012) to distinguish them from the faster counterparts.

A proposed explanation for the large spatial velocity of runaways is that they originate from disrupted binary systems (Zwicky 1957; Blaauw 1961; Boersma 1961), which naturally explains the lower number of companions they have, compared to typical massive stars (e.g., Blaauw 1961; Gies & Bolton 1986; Sana et al. 2014). Figure 1 sketches the typical evolution of a massive binary system. Most binaries are disrupted at the time of the first core-collapse.

Fig. 1. Schematic depiction of the evolution of a close massive binary through stable mass transfer. The evolution forks at the first core-collapse event in the system: the vast majority of the systems are disrupted and each produces a runaway or walkaway star which can travel far from its birth location. The systems remaining bound are possible progenitors of X-ray binaries and gravitational wave sources. The fraction of binaries disrupted 𝒟 comes from the simulations presented in this study.

An alternative mechanism to produce stars with peculiar spatial velocities is dynamical ejection from a star cluster (e.g., Poveda et al. 1967; Leonard 1991). Dynamical interactions with a supermassive BH can also disrupt a binary, but the ejection velocities achieved in this scenario are typically much higher (≳10³ km s⁻¹, Hills 1988).

Both mechanisms, i.e., the “binary disruption scenario” and the “dynamical ejection scenario”, are expected to act in nature, but their relative importance is not not well constrained (e.g., Hoogerwerf et al. 2000, 2001; Guseinov et al. 2005). Hoogerwerf et al. (2001) analyzed the properties of a sample of 56 nearby runaway stars and 9 radio pulsars, and traced back the runaways to their most likely parent stellar group. From the sub-sample for which a clear identification of the parent group was possible, they estimated that the disruption of binaries is responsible for roughly two thirds of observed runaways (see also Gies & Bolton 1986; Gies 1987; Stone 1991; Hoogerwerf et al. 2000; Dinçel et al. 2015; Boubert et al. 2017a). However, this claim could not be confirmed in the re-analysis of the same sample by Jilinski et al. (2010), who found that most of the runaways in the sample were bound spectroscopic binaries.

Unbound stars resulting from the disruption of a binary are of potential interest for several topical questions in astrophysics. Their kinematics and stellar properties bear imprints of the uncertain physical processes that govern the evolution of their binary progenitor systems, including the expected phase of mass transfer between the two stars (phase B. in Fig. 1). Detailed knowledge of the spatial velocity of runaway stars can improve the accuracy of wind mass-loss rate determinations relying on their bow shocks (e.g., Gull & Sofia 1979; Kobulnicky et al. 2018). Of particular interest is the question of whether they can provide any unique constraints on the physics of core-collapse, in particular the natal kick on the compact objects that they produce. Such kicks are expected either from asymmetries in the explosion and/or the neutrino emission (Shklovskii 1970; Wongwathanarat et al. 2013; Janka 2013, 2017), and they determine which systems disrupt and eject the companion star and which systems remain bound and thus have a potential as future X-ray and gravitational wave sources.

These unbound stars have also been considered for their potential importance as non-canonical sources of stellar feedback (e.g., Larson 1974; Ceverino & Klypin 2009). They can travel long distances and end their lives tens to hundreds of parsecs away from their birth location. Because of their motion, their ionizing photons are less likely to be absorbed by their birth clouds. Therefore, in the context of the re-ionization of the early universe, the ionizing radiation of stars ejected from a binary are more likely to escape and become available for ionization of hydrogen in the interGalactic medium (e.g., Conroy & Kratter 2012; Kimm & Cen 2014; Ma et al. 2016). Furthermore they are expected to explode in lower density regions, which may change the impact they have as sources of turbulence in the ISM (e.g., Gatto et al. 2015).

The present and upcoming Gaia data releases provide an important observational motivation for this study (e.g., Perryman et al. 2001). The determination of distances, proper motions and radial velocities by the Gaia satellite is expected to drastically increase the available sample of massive stars with precisely known velocities (Gaia Collaboration 2016a,b).

In this study, we present a systematic theoretical study of the kinematical signatures of stars ejected from massive binary systems and how these depend on the uncertain physical processes governing binary stellar evolution. For this purpose we employ a rapid binary population synthesis code, which we have updated to account for the main relevant processes that affect the disruption of the binary systems (see Sect. 2).

We first provide insight in our simulations by presenting the results for an individual system in Sect. 3. We present in Sect. 4 an analytic estimate of the typical velocity of ejected companions expected from the disruption of binaries. Then, in Sect. 5 we describe the results of a fiducial simulation of a full population of binary stars, for which we adopt realistic input assumptions provided by detailed simulations and/or observations when available. In addition, we present in Sect. 6 an extensive grid of models where we vary the most relevant uncertain assumptions. This allows us to study (i) which predictions are robust against the model uncertainties and (ii) which uncertain physical processes leave a clear imprint on the observables that can be used to test and constrain them.

Various earlier studies have discussed the evolution and interaction of populations of massive binary systems. Examples include, but are not limited to, Vanbeveren (1982), De Donder et al. (1997), Fryer & Kalogera (1997, 2001), Fryer et al. (2012), Belczynski et al. (2012), Repetto et al. (2012), Fragos et al. (2013), Grudzinska et al. (2015), Boubert & Evans (2018). Most of these studies focus on the minor fraction of systems that remain bound and are the progenitors of X-ray binaries and/or gravitational wave sources. A few studies focus specifically on disrupted systems and unbound stars. Examples are van Rensbergen et al. (1996), De Donder et al. (1997), Dray et al. (2005), Eldridge et al. (2011), Bray & Eldridge (2016), Boubert et al. (2017a), Zapartas et al. (2017a), Boubert & Evans (2018). We expand on these studies by using updated physical assumptions, and by focusing on which physical processes could be constrained using the Galactic population of unbound stars. A comparison to previously published results is given in Appendix A.

Our first main result is that the majority of disrupted binaries eject a slow-moving walkaway star (v <  30 km s⁻¹), rather than a faster runaway. This result is robust against variations of uncertain parameters in the model, as we discuss in Sect. 6. This would imply that the more easily detected runaway stars only reveal a small subset of the population of unbound stars that are former companions of disrupted binary systems. We discuss the implications in Sect. 8.

Our second main result is that the kinematic properties and absolute number of unbound companions depend sensitively on assumptions concerning the natal kicks of BHs. We show the imprints this leaves on the mass distribution and discuss whether future observations can be used to constrain these processes.

Section 7 briefly describes the population of systems remaining bound after the first core-collapse.

We also find that the fraction of runaways among massive stars predicted by our simulations is much lower than the observed value in all our parameter variations. This finding agrees with the results of Eldridge et al. (2011), but is in potential contrast with the observational results of Hoogerwerf et al. (2001), and needs further investigation.

Finally, we discuss how our results could be used to include the effect of unbound binaries in models for stellar feedback. We provide our numerical results online¹ and at the CDS.

2. Binary population synthesis calculations

We carry out population synthesis calculations of isolated binaries with the rapid binary evolution code binary_c (Izzard et al. 2004b, 2006, 2009, 2018; de Mink et al. 2013; Schneider et al. 2015). This code is based on the algorithms by Tout et al. (1997), Hurley et al. (2000, 2002), which rely on the analytic fits to the single stellar evolution models from Pols et al. (1998).

First, we compute a fiducial population using observationally favored assumptions for the free parameters that describe the initial conditions and physical assumptions. Then, we check the robustness of our results (or equivalently, to which assumptions they are most sensitive) by varying the free parameters that are most relevant for the velocity distribution of disrupted binaries. To limit computational costs, we explore variations in each parameter one-by-one while keeping the other parameters fixed to our fiducial choices, following the approach of Fragos et al. (2013), de Mink et al. (2013, 2014), Zapartas et al. (2017a,b), Belczynski et al. (2017). Effectively, binary_c treats each parameter as independent from the others, therefore this approach does not account for possible correlations between either initial distributions (e.g., the one between initial period and mass ratio suggested by Moe & Di Stefano 2017) or uncertain physical processes (e.g., core spin and natal kick amplitude). Varying two (or more) parameters simultaneously, or changing the algorithmic representation of the uncertain physical processes, might possibly result in larger variations than those we present in Sect. 6.

2.1. Initial distributions and parameters

Each binary system in our calculations is characterized by a zero age main sequence (ZAMS) mass for the primary², initial mass ratio , and initial orbital period P^ZAMS. binary_c builds a grid in this parameter space and weighs each system according to the initial distributions described below. We present our results in terms of probability per binary system. This can be converted in a probability per unit stellar mass by dividing it by the mean mass of a binary system in the population (0.42 M_⊙ for our fiducial assumptions, assuming there are no binaries for ≤ 2 M_⊙, see below).

We select N_M1 primary stars with ZAMS mass at logarithmically spaced intervals in the range 7.5 M_⊙ ≤  ≤ 100 M_⊙. We weigh each primary star with an initial mass function (IMF) with slope α = −2.3 (Kroupa 2001). In our model variations, we explore values of α of −1.9 (Schneider et al. 2018) and −3.

For each primary star of mass , we select N_q secondaries with mass , taken at regular intervals in mass ratio q between 0.1 and 1 assuming a flat distribution (slope κ = 0, e.g., Kouwenhoven et al. 2005). We consider in Sect. 6 also variations with κ = ±1.

Finally, for each pair of masses (, ) we choose N_p different periods P^ZAMS equally spaced in log₁₀(P^ZAMS/days) between 0.15 and 5.5. We weigh the birth probability of each binary system with a mass-dependent distribution for the initial orbital period: if M₁ <  15 M_⊙ we assume a flat distribution in log₁₀(P^ZAMS/days) (Öpik 1924; Kobulnicky & Fryer 2007); while for M₁ ≥ 15 M_⊙ we assume a powerlaw distribution in log₁₀(P^ZAMS/days) with exponent π = −0.55 (Sana et al. 2012). We explore values of π = 0 and π = −1 (for all values of ) in our model variations. We chose our upper-limit on the initial period, 10^5.5 days, to include wide systems in which the stars effectively evolve as single stars (de Mink & Belczynski 2015).

To limit the dimensions of the parameter space, we assume all orbits to be initially circular (i.e., eccentricity e = 0). Because we focus on the population of stars resulting from the disruption of post-interaction binaries this assumption is not critical (de Mink & Belczynski 2015): orbits are expected to circularize because of tides before or during mass transfer (Belczyński & Bulik 1999; Hurley et al. 2002, but see also Eldridge 2009 for arguments against circularization of the orbit in post-main sequence mass transfer).

In our fiducial model, we assume a canonical metallicity (i.e., mass fraction of elements heavier than helium) Z = 0.02. This value is slightly above the most recent determination for solar neighborhood Asplund et al. (2009), but it should generally describe young population of massive stars in the Milky-Way. The models from Pols et al. (1998), which serve as input for our computations, adopt the isotopic mixture of Anders & Grevesse (1989). In Sect. 6, we also consider Z = 0.0002, 0.0047, 0.008, and 0.03.

We initialize the stellar rotation rate with a mass-dependent equatorial velocity according to Hurley et al. (2000), assuming alignment between the stellar spins and the orbital angular momentum. In one of our parameter variations, we draw the initial spin velocity randomly from the distribution given by Ramírez-Agudelo et al. (2015) (R-A15, see last line in Table 1), but maintain the assumption of alignment. However, changing the initial rotation rate has almost no impact on our results because binary interaction processes overwrite it (de Mink et al. 2013).

We build a grid in the initial parameter space, and initialize a binary system in each cell of such grid. Each system is assigned a probability which can then be multiplied by the amount of time spent by the system in that cell to obtain a quantity comparable to observed number counts. See Appendix B for more details.

For each set of parameters, the total size of our model grid is N_M1 × N_q × N_p = 50 × 50 × 100 = 250 000 binary systems³. Increasing this resolution to N_M1 × N_q × N_p = 100 × 100 × 200 introduced no significant variation in our results.

2.2. Physical assumptions

We follow the evolution of each binary system until the first core-collapse (CC) event or until they merge. Stellar winds can have a dramatic effect on the evolution of single stars (e.g., Renzo et al. (2017)) and impact significantly the orbital evolution in a binary. We include wind mass loss in our single star models as implemented by de Mink et al. (2013), i.e. we use a combination of mass loss rates from Vink et al. (2000, 2001) for hot hydrogen-rich stars, Nieuwenhuijzen & de Jager (1990) for cool stars, including a luminous blue variable (LBV) enhancement as in Hurley et al. (2000), and from Hamann & Koesterke (1998) reduced by a factor of 10 for Wolf-Rayet stars. We also include the wind mass loss enhancement for fast rotating stars following Maeder & Meynet (2000).

We include the effects of tides on the spins and orbital angular momentum using the algorithm from Hurley et al. (2002) based on the calculations of Zahn (1977) and Hut (1981).

When the radius of one star exceeds the Roche radius, (calculated using Eggleton 1983 fitting formula) we use the algorithm of Claeys et al. (2014) to determine the mass transfer rate (see their Eq. (10)).

In our fiducial simulation, we assume a variable mass transfer efficiency. We limit the accretion rate of the accretor to 10 times its Kelvin-Helmholtz timescale, i.e., we assume a mass-transfer efficiency β_RLOF = β_th. Larger rates are likely to drive the accretor out of thermal equilibrium⁴ and lead to unstable mass transfer (Neo et al. 1977; Hurley et al. 2002). For a large portion of the parameter space, β_th results in a rather conservative mass transfer (Schneider et al. 2015). Since the efficiency of mass transfer is a major uncertainty in binary evolution (de Mink et al. 2007), we also consider parameter variations with β = 0,  0.5, and 1, respectively, which bracket the range of physical possibilities.

Matter that is not accreted is assumed to leave the system with the specific angular momentum of the orbit of the accretor (e.g., Soberman et al. 1997; van den Heuvel et al. 2017), h = γ_RLOFJ_orb/(M_don + M_acc), with γ_RLOF = M_don/M_acc, where J_orb is the orbital angular momentum, M_don and M_acc the donor and accretor masses, respectively. In our model variations, we also explore a scenario where the mass that is not accreted leaves the binary from the outer Lagrangian point L2 and is assumed to form a circumbinary disk (i.e., , Artymowicz & Lubow 1994). We further consider the assumptions that the mass that is not accreted leaves the system with the orbital specific angular momentum, i.e. γ = 1, which is the standard model of Belczynski et al. (2008), Dominik et al. (2012, 2013).

If, by the time the first star fills its Roche lobe, the mass ratio of the accretor to the donor is smaller than a certain threshold (M_acc/M_don <  q_crit), we assume that the system enters a common envelope phase. The threshold value q_crit is uncertain and depends on the evolutionary stage of the donor. In our fiducial simulation q_crit, A = 0.65 for a main sequence (MS) donor (de Mink et al. 2007), q_crit, B = 0.4 for a Hertzsprung gap donor (Hurley et al. 2002), and q_crit, RSG = 0.25 for core-helium burning and red supergiant donors (Claeys et al. 2014). We also consider values of q_crit, A = 0.25,  0.8, q_crit, B = 1,  0.5,  0, and q_crit, RSG = 1.0 in our model variations.

We treat common envelope evolution using the α_CEλ-formalism (Webbink 1984; Livio & Soker 1988; de Kool 1990; Hurley et al. 2002). In all variations, we use an analytic fit to the λ_g values of Dewi & Tauris (2000) for the binding energy parameter λ. Those values do not include the energy stored in thermal motions and the ionization state of the material within the envelope. It is unknown what fraction of that internal energy is useful in unbinding the envelope (see, e.g., Dewi & Tauris 2000; Ivanova et al. 2013, and references therein). We assume α_CE = 1 for our fiducial simulation, i.e., perfectly efficient use of the liberated gravitational potential energy from the orbit but without additional energy sources. We explore variations with α_CE = 0.1, i.e., inefficient use of the orbital energy to eject the common envelope, and α_CE = 10, intended as indicative of a fairly extreme case of additional energy input. We also test the combination of α_CE = 0.1 and 10 with q_crit, RSG = 1.0, to test the influence of the efficient use of energy in common envelope ejection in cases when a common envelope is more likely. For binaries that interact while neither star has a well-defined core-envelope structure, such as MS stars and helium-MS stars, we assume that the common envelope phase leads to a merger as described in de Mink et al. (2013), Schneider et al. (2015).

2.3. Natal kick

When a CC event happens in a binary system three different physical ingredients contribute to the possible disruption of the system and the ejection of the companion star.

i) The orbit is modified by the sudden change in gravitational potential, because of the mass lost through ejecta (the so-called “Blaauw kick”, Zwicky 1957; Blaauw 1961; Boersma 1961). However, this is not the dominant effect, because the typical evolution of a massive close binary system involves mass transfer through Roche lobe overflow (RLOF) before the first CC (cf. Fig. 1). RLOF removes the envelope of the initially more massive star, limiting the amount of mass that can be ejected by the SN. The Blaauw kick alone rarely unbinds the system (Huang 1963; Tutukov & Yungelson 1973; Leonard et al. 1994).

ii) The SN shock can interact with the companion star. (e.g., Wheeler et al. 1975; Liu et al. 2015; Rimoldi et al. 2016; Hirai et al. 2018). The shock can dynamically remove mass from the companion (stripping), heat the envelope of the secondary, thus enhancing its own mass loss (ablation), and deposit mass and momentum on the secondary. The injection of energy can inflate the companion and make it look redder for a duration comparable to its thermal timescale. Our treatment of accretion only changes the mass of the companion, without checking for structural readjustments of the star. However, the interaction between the ejecta and the secondary is typically a small effect (e.g., Liu et al. 2015).

iii) The natal kick of the compact object changes its kinetic energy and momentum. This kick is caused by asymmetries in the SN ejecta and/or neutrino flux at the explosion, possibly seeded by the late core and shell burning phases of stellar evolution (see e.g., Wongwathanarat et al. 2013; Holland-Ashford et al. 2017; Grefenstette et al. 2017; Katsuda et al. 2018). It is typically parametrized using a kick velocity v_k, drawn from a distribution of amplitude v_k ≡ |v_k| and directions (see below). The natal kick is the dominant reason for the disruption of binaries.

We have updated the treatment of binary disruptions by a CC event in binary_c, following Tauris & Takens (1998; hereafter TT98). To re-calculate the post-CC orbital parameters, their algorithm assumes instantaneous loss of the ejecta, because the ejecta velocity v_ej is much larger than the orbital velocity v_orb (v_orb ≲ 10² km s⁻¹ ≪ v_ej ∼ 10⁴ km s⁻¹). This algorithm considers the motion of the compact object within the pre-CC orbit, and the interaction between the ejecta and the companion star. We set the relative change of mass of the secondary due to stripping, ablation, and accretion of SN ejecta following a fit to the three-dimensional hydrodynamical simulations of Liu et al. (2015) for the impact of the ejecta on a M₂ = 3.5 M_⊙ star, which is the most massive companion they considered. We also assume an efficiency of momentum transfer⁵ from the ejecta to the companion star of 0.5 (TT98). The contribution of the interactions between the SN ejecta and the secondary star to the total natal kick is typically small, ≲10 km s⁻¹ (TT98; Liu et al. 2015; Hirai et al. 2018).

In our fiducial run, the kick direction is assumed to be isotropically distributed in the frame of the collapsing star (Wongwathanarat et al. 2013; Bear & Soker 2017), although Johnston et al. (2005) and Kaplan et al. (2008) suggest that there is weak evidence for kick-spin alignment in well observed pulsars such as the Crab or Vela. Johnston et al. (2005) also underlined that the direction of the pulsar spin does not necessarily match the direction of the pre-CC spin of the stellar core, if significant torques develop during the CC (e.g., Kazeroni et al. 2016). In our model variations, we also consider kicks constrained in a cone with an opening angle of α = 10° oriented along the spin of the exploding star (which we assume to be perpendicular to the orbital plane), and kicks constrained at angles from the orbital plane 90 − α ≤ 45°. The kick amplitude |v_k|≡v_k is drawn from a Maxwellian distribution with one-dimensional root mean squared dispersion σ_kick = 265 km s⁻¹ (Hobbs et al. 2005). Such distribution is motivated by the observation of the proper motions of pulsars in Hobbs et al. (2005) (see also Lyne & Lorimer 1994). We also compute populations with σ_kick = 0,  300,  1000 km s⁻¹ in our model variations.

The value of v_k drawn is then reduced to consider the amount of matter that falls back after the successful launch of the SN shock, i.e. v_k → v_k(1 − f_b), were f_b is the fallback fraction taken from⁶ Eq. (16) in Fryer et al. (2012), corresponding to their “rapid SN engine” which reproduces the NS-BH mass gap between ∼2 − 5 M_⊙ (e.g., Farr et al. 2011, however see also Wyrzykowski et al. 2016 regarding the existence of this mass gap). The inclusion of fallback also determines the mass of the compact remnant obtained after each CC event, and we set the mass boundary between NSs and BHs to 2.5 M_⊙.

The fallback fraction is highly uncertain in BH formation. The algorithm we use here assumes a large amount of fallback, implying close to zero BH kick amplitudes, sometimes referred to as a “BH momentum kick” (Belczynski et al. 2008; Stevenson et al. 2017). Whether this is realistic is subject of debate in the literature. Evidence for non-zero BH natal kicks comes from the observed Galactic latitude of BH X-ray binaries (e.g., Fragos et al. 2009; Repetto et al. 2012, 2017; Repetto & Nelemans 2015, although see also Mandel 2016), the possibility of retrograde BH spin (e.g., Morningstar et al. 2014), eccentric orbits in these systems (e.g., Remillard & McClintock 2006), the small number of Wolf-Rayet star-BH binaries (e.g., Dray et al. 2005) and the gravitational wave constraints on BH spins (e.g., O’Shaughnessy et al. 2017; Wysocki et al. 2018). Recently, multi-dimensional calculations of fallback in CC resulting in the formation of a BH found that large natal kicks might be possible (Chan et al. 2018). Therefore, we also consider a model without the rescaling of the kick amplitude, but still including fallback for the remnant mass calculation. This variation is effectively equivalent to a so-called “velocity kick” for the BHs. To test an intermediate BH kick, we run a simulation using σ_kick = 100 km s⁻¹ and no fallback down scaling of the kick amplitude for the CC resulting in BH formation, and our fiducial kick (σ_kick = 265 km s⁻¹, including fallback scaling) for CC resulting in NS formation.

Another debated issue is whether the CC of low-mass iron (or oxygen-neon-magnesium) cores produces small natal kicks, possibly resulting in a bi-modal kick distribution (Katz 1975; Arzoumanian et al. 2002; Pfahl et al. 2002; Podsiadlowski et al. 2004; Knigge et al. 2011; Beniamini & Piran 2016; Tauris et al. 2017; Verbunt & Cator 2017; Verbunt et al. 2017). To explore the possibility of low kicks for low mass collapsing cores, we run one model variation reducing the natal kick for NSs less massive than 1.35 M_⊙ (Schwab et al. 2010; Knigge et al. 2011) drawing the natal kicks for these NSs from a Maxwellian with σ_kick = 30 km s⁻¹ creating a double peaked distribution. For more massive NSs and BHs the kick is drawn from a Maxwellian with σ_kick = 265 km s⁻¹.

In all our simulations, we draw N_kick = 20 intrinsic natal kick directions and amplitudes for each CC in our population, effectively increasing the size of our sample by this factor. Convergence tests showed no significant variations of our results increasing N_kick to 50.

2.4. Caveats

In addition to the limitations discussed earlier in this section there are a few more caveats that should be kept in mind when comparing our simulations with observations (aside from the possible effects of biases).

We only consider runaways and walkaways resulting from the disruption of binary systems at the first stellar collapse. Observational samples of runaway stars also include another sub-population, coming from cluster ejection (e.g., Poveda et al. 1967; Leonard 1991) Similarly, walkaway stars produced by the disruption of binaries will add to slow-moving stars generated by the dissolution of clusters which did not go through binary interactions, see e.g. Allison (2012). In Sect. 8, we consider how spectroscopic measurements may allow to disentangle binary products from the outcome of purely dynamical processes.

All the calculations presented here implicitly assume a constant star formation history (SFH) for a duration longer than the longest stellar lifetime of interest. However, a recent change in the SFH of a region will affect the ratio of runaways and walkaway stars over the number of normal stars. For example, a recent increase in the SFH leads to an increase of systems that did not yet have time to interact. For a binary to be disrupted, a CC event is needed, and it can only happen after the lifetime of the primary star is over. This introduces a delay corresponding to at least the lifetime of the most massive stars (∼3 Myr, Zapartas et al. 2017b), between a SFH event and the first binary disruption. A recent increase in the SFH can therefore decrease the fraction of runaway and walkaway stars in a population.

Finally we warn about the possible effects of stochasticity (Justham & Schawinski 2012; Eldridge 2012). We have aimed to simulate sufficient binary systems such that our results are converged, but it should be kept in mind that the fastest runaways come from relatively rare channels. In a relatively low mass stellar population one will not fully sample these rare events and stochastic effects can be large.

3. Example of the evolution of a binary system

We first describe the typical evolution of a binary system that produces an unbound companion as a result of the disruption of the system by the first CC event. The aim is to provide the reader with some insight in to our simulations. As a representative example, we consider a = 20 M_⊙ primary star, a = 15 M_⊙ secondary with P^ZAMS = 100 days (corresponding to a separation a_ZAMS ≃ 300 R_⊙). Figure 2 shows the evolution of the masses, orbital period, and orbital velocity of the secondary (and ejection velocity, after the CC of the primary) as a function of time. The labels on the top axis correspond to the phases depicted in the cartoon shown in Fig. 1.

Fig. 2. Evolution of an example binary system with = 20 M⊙, = 15 M⊙ with PZAMS = 100 days, described in Sect. 3 (see also the cartoon in Fig. 1). Top panel: primary, secondary, and total mass. Bottom panel: orbital velocity of the secondary (left y-axis) and period of the binary (right y-axis). The labels indicate the following phases: A. the MS of the primary, B. case B RLOF, emphasized by the vertical dashed line, C. the binary is made of a He star (possibly looking like a Wolf-Rayet star) and a rejuvenated MS secondary; D. (and the vertical dot-dashed line) marks the explosion of the primary in a SNIb/Ic which unbinds the NS and its companion. During the phase E1. the secondary is a slow moving walkaway star, traveling for a duration of about half the main sequence lifetime of the primary.

3.1. Binary evolution until first CC

Our example system remains detached during the MS evolution of the primary (phase A. in Fig. 2). Because of wind mass loss, the period increases by about 7%.

After ∼8.8 Myr, the primary leaves the MS, and starts expanding on a thermal timescale as it burns hydrogen in a shell. Shortly after, at the vertical dashed line B. in Fig. 2, the primary fills its Roche lobe initiating case B RLOF (Kippenhahn & Weigert 1967). For a system which has initial mass ratio close to one (q_ZAMS = 0.75 in this example), RLOF is almost conservative in our calculations, i.e. nearly all the mass lost by the primary is accreted by the companion (cf. Schneider et al. 2015). The secondary star accretes ΔM₂ ≃ 14 M_⊙, and becomes the more massive star, now with a mass of about 28 M_⊙ (red line in the top panel of Fig. 2). It is also spun up to near to breakup rotation (Packet 1981; de Mink et al. 2009), and it is rejuvenated because of the increased mass and consequent growth of its convective core (Hellings 1983; Schneider et al. 2016). The primary becomes a ∼5.3 M_⊙ helium star (Götberg et al. 2017, 2018), and the orbits widens during mass transfer, reaching a period longer than 600 days (phase C. in Fig. 2).

After about 10 Myr, the primary reaches CC (dot-dashed vertical line D. in Fig. 2). Just before the core-collapse of the primary, the orbit is circular with a period of about 700 days. The secondary has a mass of about 28 M_⊙ and a pre-SN orbital velocity of about v₂ = (M₁/(M₁ + M₂))v_orb ≃ 12.5 km s⁻¹.

If the system is disrupted at the first CC event, the final spatial velocity of the secondary is nearly equal to its pre-explosion orbital velocity (e.g., Blaauw 1961; Eldridge et al. 2011). For comparison, similar initial masses but a shorter initial period of P^ZAMS = 7 days would experience a very similar evolutionary path, but it would only widen to a pre-CC period of ∼45 days, resulting in a pre-CC orbital velocity of the secondary slightly larger than 30 km s⁻¹.

Because of the rejuvenating effect of mass accretion, the secondary remains on the MS for another ∼3.6 Myr after the CC of the primary.

3.2. Effects of the natal kick on the post-CC orbit

In our example system, the CC of the primary results in a stripped-envelope SN of type⁷ IIb/Ib/Ic (Filippenko 1997), which forms a NS of about M_NS = 1.6 M_⊙, with a fallback fraction f_b = 0.13. The SN ejecta mass is , where is the primary mass just before the CC, and M_NS is the mass of the resulting NS. Because the ejecta mass is less than half of the total pre-CC mass of the system, , the loss of the SN ejecta alone is insufficient to unbind the system (Blaauw 1961). Whether the system is disrupted depends on the amplitude and direction of the natal kick imparted to the newly formed NS.

We distinguish three cases: (i) a large kick with v_k ≫ v_orb results in disruption of the binary regardless of the direction, (ii) an intermediate kick v_k ≃ v_orb disrupts the system only if the kick orientation is favorable, while (iii) a small kick, v_k ≪ v_orb, never unbinds the binary, although it can still modify its orbit (e.g., Brandt & Podsiadlowski 1995; Kalogera 1996, TT98).

To illustrate each case, we consider here three different fixed natal kick amplitudes v_k ≃ 870, 230, and 43 km s⁻¹ (respectively corresponding to 1000, 265, and 50 km s⁻¹ before the rescaling of the kick amplitude with the fallback fraction). For comparison, the pre-CC orbital velocity in our example system is km s⁻¹ (corresponding to km s⁻¹). For each kick amplitude v_k, we draw 10 000 different directions, isotropically distributed in the frame of the exploding star (e.g., Wongwathanarat et al. 2013).

For large natal kicks, our example system is disrupted in 99.9% of all the directions drawn. In the remaining 0.1%, the compact object is shot into the envelope of the secondary (see Sect. 7 for a brief discussion). In the case of disruption, the ejected star acquires a velocity km s⁻¹. With v_dis <  30 km s⁻¹, it is below the threshold to be observationally classified as runaway, and thus it is an example of what we refer to as a walkaway star.

In the case of intermediate kick amplitudes (here we use 230 km s⁻¹), the angle θ between the kick direction and the pre-CC orbital velocity of the primary, given by , determines whether the kick disrupts the system. Kicks that are roughly aligned with the pre-CC orbital velocity, i.e. 0 ≲ θ ≲ π/4, successfully disrupt the binary. Conversely, kicks that are roughly oriented opposite to the pre-CC orbital velocity of the exploding star just tend to slow the orbital motion: the system will remain bound on an eccentric orbit. In our example, we find that about 84% of the kick directions result in disruptions and produce a walkaway star with .

With small kick amplitudes (here v_k ≲ 50 km s⁻¹), our example system is never disrupted. The post-SN eccentricity, separation, and systemic velocity of the binary vary with the kick direction (e.g., Brandt & Podsiadlowski 1995; Kalogera 1996; Tauris & Takens 1998). We discuss the population of binaries surviving the first CC (and hosting a compact object) in Sect. 7.

In wide pre-CC systems, small kicks have a larger relative impact on the velocity of the unbound secondary. The reason is that in wide pre-CC orbits, smaller kicks result in a longer time for the compact object to exit the orbit and exert a gravitational pull on the other star. Conversely, for short pre-SN period systems (unlike the example considered here) the effects of large kicks are more important. This is because the time for the compact object to exit the pre-CC orbit is short, but larger kicks are more likely to result in the disruption of the system than small kicks.

4. Analytic estimates

The key quantity to define runaway and walkaway stars is their velocity. As the example described above demonstrates, the final velocity obtained at the binary disruption is of the order of its final (pre-CC) orbital velocity (Blaauw 1961; Eldridge et al. 2011). Therefore, using a few simplifying assumptions, we can estimate the ejection velocity for systems that have experienced stable mass transfer, and express it in terms of the initial parameters of the binary. In particular, we will express the final orbital velocity as a function of the initial orbital velocity in the binary, and use it as a proxy for the ejection velocity. This gives insight on the possible outcomes of the binary evolution since we can estimate the typical velocity at ZAMS from observational constraints on the binary populations, (e.g., Sana et al. 2012).

Using Kepler’s third law, we can write the orbital velocity v₂ of the secondary star in the frame of the center of mass at any time during the binary evolution as

(1)

where G is the gravitational constant, a is the semimajor axis of the orbit (assumed to be circular), and M₁, M₂ are the two masses.

If we assume that a constant fraction β_RLOF = |Ṁ_acc|/|Ṁ_don| of the mass transferred is accreted by the secondary star, and that the mass not accreted leaves the binary system with a constant fraction γ_RLOF of the total specific orbital angular momentum, we can relate the final and initial separations with

(2)

The parameter β_RLOF enters implicitly in the values of the final pre-CC masses (see Eq. (3)). Equation (2) neglects for simplicity all other processes that can remove mass from the binary, such as stellar winds (see Soberman et al. 1997, for a more elaborate expression).

The most common type of binary interaction is case B RLOF (Kippenhahn & Weigert 1967). During stable case B RLOF, the donor star typically loses nearly all its hydrogen-rich envelope (Götberg et al. 2017, 2018; Yoon et al. 2017). We parametrize the fraction of the total mass in the envelope as μ_env = M_env/M, or in other words 1 − μ_env is the fraction of mass in the helium core. The parameter μ_env is dependent on the overshooting assumed in the underlying stellar evolution models. The single massive star models from Pols et al. (1998) give values in the range 0.6 ≲ μ_env ≲ 0.8 with the larger values corresponding to less massive (and thus more common) stars. More recent models (e.g., Brott et al. 2011) typically adopt larger values for the overshooting parameter, resulting in smaller values of μ_env. Using this parametrization, we can relate the final masses of the two star in the binary to their ZAMS masses as

(3)

We can now express the final orbital velocity of the secondary in terms of the initial . Introducing , β ≡ β_RLOF, μ ≡ μ_env, and γ ≡ γ_RLOF to simplify the notation, we can write

(4)

In Fig. 3 we visualize the outcome of these analytic estimates for different assumptions for the mass transfer efficiency β ≡ β_RLOF and the angular momentum loss parameter γ ≡ γ_RLOF. To make this figure, we adopt the parameters of the example system of Sect. 3 (q = 0.75,  μ_env = 0.75).

Fig. 3. Analytic estimates of the final-to-ZAMS ratio of the specific orbital angular momentum h (left panel), the orbital separation a (center panel), and orbital velocity of the secondary v2 (right panel), assuming qZAMS = 0.75 and μenv = 0.75, which are representative for the example system of Sect. 3. Blue colors correspond to a decrease, red to an increase. The mass transfer efficiency βRLOF and the angular momentum parameter γRLOF are assumed constant throughout the evolution.

The left panel of Fig. 3 shows the relative variation in the specific orbital angular momentum . For fully conservative evolution (β = 1, white region on the right of this panel), h stays constant. Also assuming γ = 1 (e.g., Dominik et al. 2012) keeps the specific angular momentum constant, since this corresponds to losing mass and angular momentum at the same relative rate. Conversely, low values γ ≲ 1 result in a net increase of the specific orbital angular momentum because mass is lost relatively faster than angular momentum, while very large γ ≳ 1.5 values result in a loss of specific angular momentum.

The central panel of Fig. 3 shows how the evolution affects the ratio of the final-to-initial separation (cf. Eq. (2)). The vast majority of the parameter space results in significant orbital widening (red colors), with only rather extreme angular momentum losses (γ ≳ 2) in non-conservative (β ≲ 0.6) systems result in orbital shrinking.

The right panel in Fig. 3 shows the main parameter of interest, the ratio of the final-to-initial orbital velocity of the initially less massive star in the system (cf. Eq. (4)). Most assumptions for β and γ result in a decrease in the orbital velocity of the secondary (blue colors). This is the combined effect of (i) the orbital widening, and (ii) the increase mass of the secondary because of mass transfer. Only for very non-conservative systems (β ≲ 0.3) experiencing large angular momentum losses (γ ≳ 4, corresponding to non-accreted mass removing four times the specific orbital angular momentum of the binary) the final velocity of the secondary is higher than its initial value.

In our numerical simulations we adopt physically motivated model for the efficiency of mass transfer instead of constant parameters, β and γ, as we assumed in these analytical considerations. As we will argue later, most of the massive unbound companions result from nearly conservative mass transfer. Therefore, it is instructive to consider the analytical solution for β = 1. The expression for the final orbital velocity of the secondary (Eq. (4)) simplifies to

(5)

where we use as a typical value μ = 3/4 = 0.75 in the last step. Equation (5) gives a final-to-initial orbital velocity ratio for the secondary always smaller than one for initial mass ratios q >  0.25. In other words, for all systems with initial mass ratios of interest evolving through stable and nearly conservative mass transfer, the secondary always slows down substantially, by up to a factor of 5 for an initially equal mass system. We emphasize that we have ignored the additional effect of mass loss by stellar winds.

These cause the orbit to widen further and lead to a further slow down of the secondary star.

Using the typical values of our example system discussed in Sect. 3 (q = 0.75), we find km s⁻¹. Our numerical model for this system gives an even lower value, primarily because of the additional effect of widening as a result of mass loss through stellar winds. As a sanity check we recomputed the evolution of the example system, but artificially switching off the stellar wind mass loss for both stars. This gives km s⁻¹, in good agreement with the analitycal estimate. The remaining difference is due to the fact that the mass transfer is not fully conservative.

For the analytic estimates presented in this section, we only considered the case of stable mass transfer. Binary evolution through unstable mass transfer takes place for a limited range of the initial distributions (e.g., Soberman et al. 1997; Schneider et al. 2015). It is expected to result in a common-envelope phase, followed by substantial shrinking of the orbit either leading to a merger or the formation of a compact binary system if the envelope is ejected successfully. While this is a very interesting pathway to create very fast runaway stars, the numerical simulations presented later indicate that this channel produces primarily low mass runaway stars. This is because it typically concerns system with extreme initial mass ratios in which the secondary does not significantly gains mass. This means that the secondaries in these systems are typically not very massive and therefore this channel does not significantly contribute to the production of massive early type runaways, at least in our simulations. We return to this in Sect. 8.

5. Ejected companions

In this section, we present the results obtained from our numerical simulations for a full population computed with our fiducial assumptions (cf. Sect. 2). We describe the robustness of our findings against model assumptions in Sect. 6. Our main goal is to characterize the velocity acquired by the companion stars ejected at the time of the first CC in the binary. Therefore, we ignore all binaries that merge, which constitute ∼22% of our fiducial population (cf. Fig. 4), in good agreement with Sana et al. (2012), de Mink et al. (2014), Zapartas et al. (2017a,b); see also Kochanek et al. (2014).

Fig. 4. Overview of the binary evolution scenarios up to the first CC event. The branching ratios shown are from our fiducial simulation, and we highlight in red the disruption fraction 𝒟. The errors on each fraction exclude the run without SN kicks (σkick = 0 km s−1), which produces an unrealistically low disruption fraction (cf. Table 1 and Sect. 6).

Once the merging systems are excluded, the first CC happens in the presence of a companion. The CC can result in either the disruption of the binary (𝒟 = 86% in our fiducial run), or the newly formed compact object remains bound to the companion star. In Sect. 5.1 we focus on the population of ejected MS companions, which are 75% of all the ejected companions in our fiducial run. We consider companions that have already evolved off the MS at the time of the first CC in Sect. 5.3.

For completeness, we describe the systems with a MS companion remaining bound in Sect. 7.

5.1. Ejected main sequence companions

The fraction of CC events occurring with a MS companion in our fiducial run is ∼76% (and it ranges from 69 − 90% in our parameter variations). The vast majority of these binaries which host a MS companion are disrupted at the time of the first CC. Figure 5 shows the distribution of velocities of the ejected companions. The three different colors subdivide the population based on the mass of the ejected companion. The gray histogram shows all ejected MS companions. The blue histogram shows only massive MS companions, , where is the mass of the secondary just after the collapse of the primary. This group roughly corresponds to all secondaries that will experience CC at the end of their lifetime. The red histogram shows the velocity distribution of companions with mass larger than , corresponding roughly to O-type stars. In the rest of this section, we quote the ratios for the population of massive stars (i.e. ) unless otherwise stated, see also Table 1 for more information.

Fig. 5. Velocity distribution of MS stars ejected from a binary system at the time of the first CC. Top panel: corresponding cumulative distributions. The walkaways are one order of magnitude more numerous than the runaways. The gray, blue, and red histograms show all MS secondaries, only MS stars more massive than 7.5 M⊙ (i.e., roughly those that might experience CC), and only MS stars more massive than 15 M⊙ (i.e., roughly the O-type stars), respectively. The mass considered here is taken right after the CC of the primary star. The dashed lines show the distribution of walkaways and runaways that have gone through RLOF (or common envelope evolution) before being ejected. Almost all massive runaways and walkaways, and the majority of the ejected stars at large velocities, have gone through RLOF. The inset plot magnifies the runaway regime vdis ≥ 30 km s−1. See also Fig. C.2 for a wider velocity range, Fig. C.1 for a figure accounting for the remaining lifetime as walkaway or runaway star, and Fig. 12 for an example of how the velocity dispersion of a star forming region affects these distributions.

Although we do produce runaways, the majority of the systems eject a star slower than 30 km s⁻¹, i.e., a walkaway star. We find that the CC of a star with a MS binary companion is about twenty times more likely to produce a walkaway than a runaway. In the first line of Table 1, we list the ratio ℛ of walkaway to runaway stars produced per CC event in a binary. From our fiducial simulation, we obtain ℛ_7.5 ≃ 20 for massive companions. As shown in Table 1, considering progressively larger masses, the number of walkaways produced per each runaway from a binary increases. In Sect. 8.3, we address the impact of the velocity dispersion within a star forming region on this result.

Runaways resulting from the disruption of binaries rarely exceed 60 km s⁻¹: about 99.8% of the massive unbound companions we simulate are slower than this threshold. We note that higher velocity are not strictly forbidden, but extremely unlikely (cf. Fig. C.2 which shows the velocity distribution on a logarithmic scale). Low mass runaways can reach much tighter pre-CC orbits through common envelope evolution (without significant accretion), and thus reach higher ejection velocities, but typically not in excess of 400 km s⁻¹. Massive stars significantly faster than this are likely produced by different mechanisms (see Boubert et al. 2017a; Brown 2015 for a review). However, in extremely rare cases where a very short pre-CC binary is formed and also successfully disrupted, much higher velocities can be achieved. The absolute maximum velocity we obtain is about 1110 km s⁻¹, in good agreement with Tauris (2015).

The slowest stars come from the widest pre-CC binaries: the low-velocity drop-off of the distribution in Fig. 5 is thus an effect of the upper-limit on our period distribution (corresponding to an upper-limit on the pre-CC period).

The velocity distribution peaks at v_dis ≃ 6 km s⁻¹. Ejected MS companions with such velocity typically originate from systems with q_ZAMS ≳ 0.5 and long initial periods (log₁₀(P^ZAMS/days)≃3), regardless of the stellar masses. These systems experience late case B or case C mass transfer (Lauterborn 1970), but owing to the sufficiently high mass ratio, the mass transfer remains conservative and the system widens.

Among low and intermediate mass walkaways with v_dis ≃ 6 km s⁻¹, there is also a contribution from wide systems which do not experience mass transfer. This contribution is almost negligible for ejected stars more massive than 15 M_⊙. This can be seen from the dashed distributions in Fig. 5, which show that the majority of massive runaways and walkaways have accreted mass from their companion before its CC. Over the entire velocity range shown (i.e., considering both runaways and walkaways), 71% of ejected MS secondaries more massive than 7.5 M_⊙, and 91% of those more massive than 15 M_⊙ have accreted from their companion in the previous evolution. Among low mass ejected companions, 50% have experienced either RLOF or common envelope evolution before their ejection from the binary. From the inset in Fig. 5, where the blue and red dashed lines representing post-interaction systems overlap with the full distribution, it is clear that all massive runaways from binary disruption have experienced mass transfer.

The majority of the ejected massive MS companions are not usual single stars, and this remains true for a very large fraction of the low mass MS ejected companions. Before the binary disruption, they have accreted mass from their companion. Mass transfer also causes the convective core of the accretor to grow, resulting in rejuvenation of the accreting star (e.g., Hellings 1983; Schneider et al. 2016). The fact that most ejected MS stars are accretors is the result of a combination of assumptions. For disruption to be possible, the primary needs to be massive enough to collapse. We assume a flat mass ratio distribution, implying that on average , and we also assume that mass transfer is unstable and leads to mergers if the mass ratio is too extreme (M_acc/M_don <  q_crit), therefore binaries that can eject a star have initial mass ratios closer to one than an average system. Finally, when stable mass transfer occurs, our fiducial assumption for the accretion efficiency results in rather conservative mass transfer over most of the period range considered (Schneider et al. 2015).

To make an estimate of the ratio of walkaways per runaway existing at a given time requires us to consider the finite MS lifetime of the ejected stars. Assuming that the star formation rate in the Galaxy is constant and therefore that the Galactic population is in a steady state, the ratio of walkaways per runaway is ∼13 for masses larger 7.5 M_⊙.

Accounting for the durations of each evolutionary phase, we can also quantify the fraction of MS stars more massive than 15 M_⊙ that are runaway or walkaway stars. In our fiducial simulation, we find a runaway fraction for masses larger than 15 M_⊙ of about ∼0.5% and a walkaway fraction of ∼10%. We only simulate binaries with primaries more massive than 7.5 M_⊙, therefore we do not have in our population intermediate-mass binaries which could contribute through mass transfer and/or mergers to the normalization factor for the runaway and walkaway fraction at lower masses (Zapartas et al. 2017b). We show in Fig. C.1 the distribution in velocity accounting for the lifetime of the stars in each bin. This can be directly compared to observations assuming that the effects of the Galactic potential can be neglected.

Figure 6 shows the distribution of distances traveled by massive, unbound stars, calculated by multiplying their velocity v_dis by their remaining lifetime on the MS, which includes rejuvenation if the star has accreted mass, plus 10% of the MS duration of a star with the same helium core mass, to account for the duration of the helium core burning phase in the rejuvenated star⁸. This approach neglects the effect of an external potential on the trajectory of the stars (see Boubert et al. 2017b; Boubert & Evans 2018, for how this could be done), so it is effectively an upper limit to the distance they can travel. Considering both runaway and walkaway stars with masses M_dis ≥ 7.5 M_⊙, we find the mean distance ⟨L⟩ they travel before experiencing CC to be 126 pc, see the blue vertical dashed line in Fig. 6 and also Table 1. This number rises to ⟨L_run⟩=584 pc considering only runaways (faster than 30 km s⁻¹) more massive than 7.5 M_⊙, but it remains ⟨L_walk⟩=103 pc considering massive walkaways only (slower than 30 km s⁻¹). For comparison, the typical size of OB-associations is on the order of tens of pc, while the Galactic thin and thick disks have a vertical scale height of ∼300 and ∼1500 pc, respectively (e.g., Gilmore & Reid 1983).

Fig. 6. Distribution of the maximum distance L that ejected MS secondaries can reach, neglecting any gravitational potential. The blue and red histograms are for ejected stars more massive than 7.5 M⊙ and 15 M⊙, respectively. The vertical dashed lines mark the mean values of the distributions.

5.2. Spin of the main sequence ejected companions

Figure 7 shows the equatorial rotational velocity for the ejected MS companions with M_dis ≥ 7.5 M_⊙. The top panels show the normalized one-dimensional distribution of ejection velocity for the whole sample (the same as the blue distribution in Fig. 5), and for the runaways only, respectively. The side panels show the normalized equatorial rotational velocity distribution for the whole sample, and for the runaways only, respectively. We do not include in Fig. 7 any projection effects due to the inclination of the stellar spin axis with respect to the line of sight, nor the wind spin down due to the evolution after the binary disruption. The combination of these two effects would systematically decrease the rotational velocities, allowing any physically possible projected equatorial rotational velocity to be observable.

Fig. 7. Equatorial velocity at the time of ejection for massive (Mdis ≥ 7.5 M⊙) companions as a function of the ejection velocity. Brighter colors indicate location of the parameter space more populated by our fiducial simulation. All the massive runaway stars spin with an equatorial velocity of ∼500 km s−1, close to breakup rotation, since they have accreted mass from their companion before being ejected. The spread is due to wind spin down before the binary disruption. We do not include the effect of projecting on the line of sight, or post-ejection wind spin down in this plot. The vertical dot-dashed line marks the threshold to define runaway stars.

The main feature at the top of the central panel (also visible in the rightmost panel) indicates that almost all the massive runaways are spinning nearly at breakup rotation at the time of the ejection, because of the relatively recent spin-up during the mass transfer phase(s).

The walkaways (v_dis <  30 km s⁻¹) instead can have a broader range of rotational velocities: this is because a subset of them comes from wide, non-interacting binaries. The right panel including both the walkaway and runaway population shows two peaks for and 500 km s⁻¹. The former is due to the combination of the accretion-induced spin up and the wind spin down before the ejection, which can be more efficient for walkaways experiencing early mass transfer and with longer delay between the mass transfer and the disruption of the binary. The former peak is not related to tidal locking, which is efficient only for tight orbits that would result in much larger ejection velocities v_dis. The latter is populated by secondaries ejected after a mass transfer phase (analogous to the peak in the equatorial velocity distribution for the runaways).

5.3. Ejected post-main sequence companions

Figure 8 shows the velocity distribution of companion stars which have evolved off the MS by the time of the first CC in the binary. All together, these correspond roughly to 25% of the population of disrupted binaries (cf. Fig. 4). Other kind of systems and evolutionary channels not considered here can also produce post-MS ejected companions (e.g., Justham et al. 2009; Zapartas et al. 2017b). We distinguish three sub-populations: hydrogen-rich non-degenerate secondaries (mostly of red-supergiant secondaries, corresponding to ∼4% of the entire simulated population, including mergers and non-disrupted binaries), white dwarfs (WD, ∼8% of the total), and helium stars which have lost their hydrogen-rich envelope because of previous binary interactions (∼4% of the total).

Fig. 8. Velocity distribution of companions ejected after their MS. Helium stars and white dwarfs (WD) can have much shorter pre-CC orbits because of their smaller radii, and thus they can reach much higher ejection velocities. Conversely, the number of non-degenerate, hydrogen-rich post-MS stars is dominated by red supergiants with very large radii, which can only exist at large separation from the exploding companion, and thus typically have much slower velocities. Top panel: corresponding cumulative distributions.

For the first group, characterized by large stellar radii, the velocity distribution resembles closely that of the MS ejected companions, with a large number of walkaways compared to the runaways (roughly corresponding to 96% of the hydrogen-rich non-degenerate ejected companions).

Ejected WDs and helium stars come from different evolutionary paths. Mass transfer during the primary MS can reduce the helium core mass thus slowing its evolution, while accretion of mass speeds up the evolution of the companion. In ∼23% of our full population (excluding stellar mergers), this leads to a reversal of the CC order, with the secondary exploding first (Pols 1994). In these cases the star that is ejected is the initially more massive, which has become, by the time of the CC, a helium star or a WD.

WDs and helium stars are what is left of the initial donor and they are thus less massive than their companion after the mass transfer phase, therefore they can reach velocities much larger than (super-)giant stars. The ejection velocities of the WDs span a range 15 ≲ v_dis/(km s⁻¹) ≲ 90, and helium stars are even faster with 70 ≲ v_dis/(km s⁻¹) ≲ 120. This difference arises because it takes longer for a star to become a WD than its helium core burning duration. WDs are ejected by systems where the rejuvenation of the secondary is not extreme. This suggests that these systems have a less conservative mass transfer phase than binaries ejecting a helium star, which corresponds to less rejuvenation of the secondary and more orbital widening, and thus lower velocity of the WD at ejection. Moreover, if the companion star at pre-CC stage is a WD, there is a longer time for the binary to widen, e.g., because of winds, than if it is a helium star.

Conversely, if the ejected star is a non-degenerate helium star, then the entire evolution of the rejuvenated secondary needs to be faster than the remaining lifetime of the primary (which now has a reduced mass): the time delay between the RLOF phase and the first CC ejecting the helium star is much shorter and there is not enough time for significant orbital widening.

Figure 8 also shows a small population of slow (v_dis <  30 km s⁻¹) helium stars. These are very massive (M₂ ≳ 30 M_⊙ after mass accretion) stars characterized by large Wolf-Rayet wind mass loss rates, which cause significant orbital widening before the first CC, resulting in the slow ejection velocity. We emphasize that Fig. 8 shows the velocity of stars ejected while in the corresponding evolutionary stage, and stars ejected during their MS will then evolve into post-MS hydrogen-rich, or possibly Wolf-Rayet stars, creating a more complicated velocity distribution of massive Wolf-Rayet stars (see, e.g., Dray et al. 2005; Eldridge et al. 2011).

6. Impact of model variations

Our predictions depend on a set of parametrized assumptions. We perform a systematic study of the impact of these uncertainties by varying the values of the free parameters one-by-one, as we describe in Sect. 2. In each variation, all the other parameters are set to their fiducial value.

Table 1 summarizes our results for binaries with a collapsing star and a MS companion. If our extreme assumptions do not change significantly our synthetic population, then the outcome is robust and independent of the particular physical process. Conversely, strong dependence on any one parameter indicates that this particular parameter could be constrained by comparing to observed samples, and therefore allow for physical tests of the process it represents. Overall, our results are robust and relatively insensitive to uncertain physical processes, at least within the framework of parametrizations we assume.

6.1. Robust predictions

The fraction of binaries that are disrupted by the first CC is always larger than 75%, with two notable exceptions. These concern the natal kick distribution assumed. With a double Maxwellian kick distribution with σ_kick = 30 km s⁻¹ for NSs less massive than 1.35 M_⊙ (and our fiducial kick for more massive remnants), the disruption fraction decreases to 𝒟 = 65% even though this changes the kick only for a small range of the parameter space. Low mass NSs are the typical result of progenitors with small total and core masses, which are favored by the IMF.

If instead we make the extreme assumption of zero natal kick in all CC (σ_kick = 0 km s⁻¹), then only 𝒟 = 16% of the binaries are disrupted. In these cases the disruption is only due to the mass loss at the moment of the explosion (Blaauw kick). Since the ejecta mass must be larger than half of the total mass of the system to disrupt the binary (Blaauw 1961), and since after mass transfer the primary is typically less massive than the companion, only very wide binaries which did not experience mass transfer might be disrupted in such a way. Therefore, in this extreme assumption, the ejected stars are on average less massive and slower compared to calculations including a natal kick. We emphasize that this parameter variation is in contrast with the observed distribution of pulsar velocities, and we include it for insight only. This variation does not produce any massive runaways, which is also in contrast with observations.

We obtain the largest disruption fraction when we do not rescale the natal kick with the fallback fraction: in this model, up to 97% of the systems are disrupted. This assumption is equivalent to assuming a “BH velocity kick” (e.g., Belczynski et al. 2008; Fryer et al. 2012): BHs and NSs receive the same kick amplitude at birth, despite the differences in their masses. Therefore, many CC events resulting in the formation of a massive BH with large fallback fractions (up to f_b = 1) still result in the disruption of the binary.

Despite the overall large disruption fraction in all our parameter variations, the predicted runaway fraction for masses larger than 15 M_⊙ is low for all our parameter variations and never exceeds a few percent. Similarly, the walkaway fraction is also robust, and is generally around 10%. In other words, the walkaways outnumber the runaways by a factor of 10–30 (smaller for larger masses), regardless of the physical assumptions we make. The most notable exceptions are computations at lower metallicity (see Sect. 6.3). Another case producing slightly smaller ℛ is fully non-conservative mass transfer (β_RLOF = 0): more systems are in a very close pre-CC orbits and the average velocity of an ejected secondary is higher than in our fiducial assumptions.

It is worth noting that variations in the common envelope efficiency (α_CE) do not have a large influence on this result. The most likely reason for this is the combination of our assumptions on the stability of late case B and case C RLOF and whether massive stars develop deep convective envelopes. The models from Pols et al. (1998), underlying the fitting formulae of binary_c, barely ascend the Hayashi track for M_ZAMS ≳ 15 M_⊙, and the q_crit, RSG assumed in our fiducial simulation is rather conservative in assessing the stability of late RLOF. Both these assumptions contribute in making the parameter space leading to common envelope, significant orbital shrinking, and subsequent ejection of a fast and massive⁹ runaway relatively small. Allowing for less stable late RLOF (higher q_crit, RSG) increases the average speed of ejected companions, and decreases the relative number of walkaways. however, extreme variations in this parameter (in combination with the common envelope efficiency α_CE) are not sufficient to reconcile our predictions with commonly-accepted inferences from observations (see Sect. 8.2).

The mean value of the mass of the walkaways, runaways, and all ejected secondaries (, , and ⟨M₂⟩, respectively) exceeds 7.5 M_⊙ in almost all our parameter variations. With fully non-conservative mass transfer (β_RLOF = 0), the secondary accretes no mass. This makes them on average less massive than in our fiducial assumptions at the moment of the first CC, and also allows for less orbital widening (the donor has to lose more mass before the mass ratio inverts), resulting in a faster population of disrupted stars with lower masses compared to our fiducial case. If we assume that case B mass transfer is unstable whenever the accretor is less massive than the donor (q_crit, B = 1), or if we take a negative slope of the mass ratio distribution (κ = −1), then the fraction of binary systems that merge after the end of the primary MS increases compared to our fiducial case. When a binary system merges, it cannot eject a walkaway or runaway. Finally, assuming no natal kicks (σ_kick = 0), the only systems disrupted are those where the SN ejecta consist of more than half of the total mass of the binary at the pre-explosion stage. This requirement naturally biases the population of ejected companions towards low masses. If we assume that the kick velocity is not scaled down by the amount of fallback (“velocity kick” for BHs), CC events from more massive primaries can also disrupt the system. On average, this produces more massive ejected secondaries compared to our fiducial case. Since the remaining lifetimes are shorter, the distance traveled by the population of ejected secondaries is also decreased, despite their higher velocity.

Our assumptions about the stability of case B RLOF (assumed to be stable if the accretor-to-donor mass ratio is larger than q_crit, B) impact the velocity and mass distribution of the ejected stars. The sequence from more stable to more unstable q_crit, B = 1,  0.5,  0.4 (fiducial value), 0 shows progressively higher ℛ indicating relatively fewer runaways, but higher average velocity ⟨v⟩ (and consequently also a larger ⟨L⟩). The increase in velocity and decrease in average mass of the ejected companions for higher q_crit, B indicates an increased importance of evolutionary channels involving a common-envelope evolution leading to tighter pre-CC orbits without significant accretion of mass from the companion.

The average distance traveled by walkaways and runaways is also a robust prediction. The entire population of ejected stars typically travels slightly farther than ⟨L⟩≃120 pc, which corresponds to an average distance for the runaways of ⟨L⟩_run ≳ 500 pc, while the walkaways alone reach ⟨L⟩_run ≳ 100 pc.

These distances are significantly shortened only assuming fully conservative mass transfer (β_RLOF = 1), or assuming higher angular momentum loss (e.g., through a circumbinary disk, γ_RLOF = γ_disk). The former results in higher masses of the ejected secondaries and thus smaller velocities and shorter lifetimes, while the latter results in more orbital shrinking and mergers. Another variation resulting in small ⟨L⟩_walk ≃ 66 pc is that with zero natal kicks (σ_kick = 0 km s⁻¹). As mentioned above, only the widest pre-CC binary can be disrupted in this case, so producing a much slower population of ejected MS stars.

The predicted rotational velocity of the ejected companions is also a prediction robust against parameter variations. The vast majority of massive (M_dis ≥ 7.5 M_⊙) main sequence runaways accrete mass from their companions, and thus spin up to critical rotation, before being ejected. Visual inspection of plots like the one shown in Fig. 7 for our fiducial simulation show no variation when changing the initial rotation rate of the stars, confirming that mass transfer in binaries overwrites the initial rotation of the accretor.

6.2. The mass function of massive runaways as probes for black hole kicks

We find that the mass distribution of runaways carries information about the BH kicks. With a sufficiently large sample, it might be used to discriminate between BH velocity or momentum kicks, at least in a statistical sense. We focus on the mass function of massive (M_dis ≥ 7.5 M_⊙) and fast runaways (v_dis ≥ 30 km s⁻¹) only, i.e., not including the more common walkaways, since we expect upcoming observational distributions to provide cleaner samples for these. Ejected companions with high masses are more likely to come from binaries with an initially very massive primary, so resulting in BH formation at the end of its evolution. The constraints on the amplitude of BH kicks from unbound companions are complimentary to those from binaries surviving the first CC, and Gaia might be able to probe both populations of bound and ejected companions (Breivik et al. 2017).

Figure 9 shows the mass distribution of runaways (v_dis >  30 km s⁻¹) more massive than 7.5 M_⊙ for three BH formation scenario. The bottom panel shows our fiducial simulation, in which natal kicks are damped by fallback (effectively, this is equivalent to a BH “momentum kick”). In this case, the 90th percentile of the massive runaway mass function is ∼25 M_⊙. In this model, the amount of fallback and the amplitude of the BH natal kicks are degenerate parameters. Because of the large fallback fraction prescribed by the Fryer et al. (2012) algorithm, the runaway mass function produced by this model is practically indistinguishable from that with no BH natal kicks: almost all the runaways produced come from NS-forming supernovae.

Fig. 9. Mass function of massive (Mdis ≥ 7.5 M⊙), runaway (faster than ≥30 km s−1) stars depends on the assumptions for the BH kick. Bottom panel: our fiducial simulation with σkick = 265 km s−1, including the fallback downscaling of the kick both for NSs and BHs. The fallback is an important effect only for BHs. Central panel: mass function for intermediate BH kicks (σkick = 100 km s−1, no fallback downscaling) and our fiducial kick for NSs. Top panel: distribution with our fiducial kick amplitude, but without fallback downscaling (i.e., BH velocity kick). The dashed lines indicate the 90th percentiles of the mass distributions.

The central panel of Fig. 9 shows the mass distribution assuming no fallback downscaling, but a smaller kick amplitude distribution characterized by σ_kick = 100 km s⁻¹ for the BHs (while we keep the same fiducial assumptions for NSs). Compared to our fiducial simulation, BHs receive larger kicks, and thus the runaways produced are also generally more massive, and the 90th percentile of the mass distribution shifts to ∼30 M_⊙.

The histogram in the top panel of Fig. 9 instead corresponds to a scenario in which the natal kick is not decreased by the amount of fallback (e.g., if the natal kick is entirely due to asymmetries in the neutrino emission). In this case BHs and NSs receive the same kick amplitude (BH “velocity kick”) and this results in a larger number of massive runaways characterized by a mass function skewed towards larger values. In this scenario, the 90th percentile of the mass distribution is as high as ∼50 M_⊙.

In summary, the tail of the mass function of massive runaways carries information about the BH natal kicks: finding runaways of large masses would support the possibility of large BH natal kicks. We do not expect large contamination of this tail from dynamical ejection of runaway stars, since dynamical channels tend to eject the least massive star among those interacting (although see also Fujii & Portegies Zwart 2011; Banerjee et al. 2012; Oh & Kroupa 2016).

Our result allow us to check also for a correlation between the mass of the ejected companions and its velocity. Figure 10 shows this comparison for the three BH kick variations of Fig. 9. As expected, there is a correlation between the maximum velocity that can be achieved because of the binary disruption, and the mass of the ejected MS companion. This can be seen by the shape of the colors in each panel. More massive companions have more inertia, or, in other words, more massive companion have slower pre-CC orbital velocity because of the mass dependent factor in Eq. (1).

Fig. 10. Mass-velocity correlation for the three BH kick variations discussed in Sect. 6.2. The brighter colors indicate where the probability per system is higher on this plane, the colors have a logarithmic scale. Left panel: our fiducial distribution, central panel: case for intermediate BH kicks (σkick = 100 km s−1 and no fallback), and right panel: large BH kicks (same as for NSs). The NS kick is effectively the same in all panels, since the amount of fallback is very small for SNe producing a NS.

From left to right, the BH kick amplitude increases (while effectively keeping the same NS kick, since the fallback fraction is small for CCSNe resulting in NS formation), and so does the maximum mass of the ejected companions. There is an sharp drop in the number of ejected companions for M_dis ≳ 70 M_⊙. There is a set of plausible reasons to explain this. First of all, stars in this mass range are intrinsically rare because of the IMF. Also, the lifetimes of very massive stars become very similar (∼3 Myr) above ∼50 M_⊙. This means that by the time the primary collapses, the secondary is also an evolved star that has finished burning hydrogen in its core, and it does not appear in the panels of Fig. 10. Finally, these stars tend to have large wind mass loss rates, which can be enhanced significantly when the star is spun up by the accretion: the combination of the wind enhancement with already large wind mass loss rates limits the growth in mass of the companions.

In an observed population, the short post-ejection lifetime would also make very massive ejected companions rare. However, the upperlimit of ∼70 M_⊙ is not strict, and it is possible to produce walkaways in this mass range through binary interactions, although at a negligible rate compared to less massive walkaways.

6.3. Effects of metallicity

We find that, other than the natal kick distribution, metallicity is the most important parameter influencing our results. Decreasing the metallicity from super-solar, Z = 0.03, to Z = 0.02 ≃ Z_⊙ (fiducial run), to Z = 0.008, Z = 0.0047 to Z = 0.0002 ≃ 10⁻² Z_⊙ a clear trend emerges from Table 1. At lower metallicity, the fraction 𝒟 of binaries disrupted by the first CC is lower (77% for the lowest metallicity, cf. 88% at the highest), but the stars ejected are on average faster. This can be seen from the average velocity ⟨v⟩, which goes from 11.2 km s⁻¹ to 21.6 km s⁻¹, but also from the progressive decrease of the ratio of walkaways per runaway ℛ, regardless of the lower mass cut. Nevertheless, walkaways outnumber runaways at all the metallicities we explore.

In our lowest metallicity run, the runaways (v_dis ≥ 30 km s⁻¹) travel, on average, as far as 706 pc from their birth location, and the walkaways reach, on average, 163 pc. The mean distance travelled by population of disrupted stars is 279 pc for the lowest metallicity.

The reduced likelihood of disrupting a system and increased velocity of the stars ejected can be understood in terms of the effect of metallicity on stellar radii. At lower Z stars of a given mass and evolutionary stage have smaller radii. This means that a given binary will enter into a mass transfer phase later, resulting in less conservative mass transfer and less orbital widening, making the binary harder to disrupt and the pre-CC orbital velocity of the secondary faster. This is corroborated by the fact that the average mass of walkaways and runaways decreases slightly with metallicity, which is indicative of a less conservative mass transfer. Moreover, at lower metallicity, the post-interaction orbital widening due to stellar winds is less important.

Another important physical effect, presently not included in our simulations, is that, at the end of the RLOF, the primary remains more massive at lower metallicity (owing again to the smaller radii, e.g., Götberg et al. 2017; Yoon et al. 2017). This should contribute to the reduction of orbital widening for interacting binaries, and make the collapsing star more massive and thus the disruption harder.

7. Systemic velocities of NS and BH binaries

Our simulations also provide predictions for the systemic velocities of binary systems that remain bound after the first CC (e.g., van Oijen 1989). The results of our simulations, which will be made available upon publication, also provide pre-CC and post-CC separations, eccentricities, and mass ratios.

The systems with a MS star and a compact object remaining bound are only of the total (see Sect. 6 for a discussion of model variations). These are of large interest since they can give rise to X-ray binaries if the compact object starts to accrete from the companion, and later possibly double compact objects. They typically originate from systems where the newly formed compact object received a small natal kick. This bias limits the extent to which these special systems can be used to infer conclusions about the properties of the general population of NS and BH. For comparison, the average effective natal kick for the compact objects formed in systems that remain bound is ∼66 km s⁻¹, much lower than the average of ∼330 km s⁻¹ (as drawn from our standard Maxwellian after rescaling to account for fallback).

Systems that remain bound also obtain a systemic velocity, because of mass lost during the supernova explosion while orbiting around the center of mass. We define the newly obtained systemic velocity v_sys as the velocity of the new center of mass (of the compact object plus the secondary star) in the frame of the old center of mass (e.g., Brandt & Podsiadlowski 1995; Kalogera 1996, TT98).

Figure 11 shows the normalized distribution of systemic velocities v_sys for different natal kick assumptions. The left (right) panel shows NS+MS (BH+MS) binaries, with the cumulative distributions in the top panels, respectively. We find that the distribution of systemic velocities are very insensitive to the model variations that we considered, except for the natal kicks. We only show the distributions for variations that show large deviations from the fiducial simulation.

Fig. 11. Normalized distribution of systemic velocities of binaries with a MS star remaining bound after the CC of the first star. The distributions in both panels are normalized to unity. Left panel: distribution for NS companions. Right panel: only binaries with a BH companion. The blue histogram represents our fiducial simulation. The orange solid line gives the distribution when the kick is not rescaled because of fallback, which changes the BH kicks, but not the NS kicks. The red distribution corresponds to a double Maxwellian kick distribution, which changes the NS kicks but not the BH kicks. Top panels: corresponding cumulative distributions.

For NS+MS binaries our fiducial simulation shows a median systemic velocity of ∼20 km s⁻¹, for nearly all assumptions that we consider. The strongest deviation is found with a run in which we adopt a double Maxwellian distribution, where low mass progenitors are assumed to form neutrons stars with very low kick velocities. This substantially increases the fraction of systems that remain bound to 1 − 𝒟 ≃ 35%. It primarily adds systems with low systemic velocities. The median of the normalized distribution for this simulation lies near 7 km s⁻¹.

In our fiducial simulations for BH+MS binaries we obtain a bimodal distribution with a large peak at 0 km s⁻¹ accounting for nearly 80% of all systems and a second broader component peaking at ∼10 km s⁻¹. The peak near zero results from progenitors that experience complete fallback (f_b = 1 according to the algorithm from Fryer et al. 2012). They receive no kick at all, and no mass is lost from the system (we neglect here, as elsewhere, the mass possibly lost to neutrinos during the CC, which is an ∼10% effect). For these systems we also find no changes in the separation and/or eccentricity after the CC of the first star. The sharp peak at v_sys ≃ 0 km s⁻¹ would likely be smeared out by the velocity dispersion of the region where the system formed.

Nearly all variations considered give indistinguishable results. However, we find large deviations from the simulation where we used the extreme assumption that compact objects get kicks of a similar velocity amplitude regardless of their mass (no fallback downscaling, orange in Fig. 11). This effectively results in much larger kicks for black holes and is sometimes referred to as “velocity kick”. Only 1 − 𝒟 ≃ 3% of systems remain bound in this simulation. The few resulting BH binaries typically have large systemic velocities: we find a broad distribution with a median at 50 km s⁻¹.

Our simulations also provide predictions for the distributions of further properties of the bound systems. If we exclude systems with f_b = 1, i.e. if we consider only NSs and BHs formed via fallback but still producing ejecta, we find an anti-correlation between the maximum post-CC separation and systemic velocity: the wider the system, the slower the maximum systemic velocity it can reach (see also Brandt & Podsiadlowski 1995).

Most systems remaining bound come from initially relatively tight pre-CC orbits. Typically, they have a reversed mass ratio (i.e., at the time of the first CC ) because of the previous mass transfer phase. In our fiducial run, for the systems remaining bound the average pre-CC mass ratio is . The majority of the bound systems, about 60% for all variations with non-zero natal kicks, have separations of less than about 500 R_⊙. This suggests that these systems will evolve through a phase during which they might be detectable as X-ray sources, and possibly even through a mass transfer phase. The tail of the distribution at large separations extends to post-CC separation as large as ∼6000 R_⊙.

We do not consider cases in which the compact object is shot within the Roche lobe of the secondary. A prompt collision between the newly formed compact object and the companion star could lead to the formation of a Thorne-Zytkow object (if the compact object is a NS, Thorne & Zytkow 1975, 1977; Leonard et al. 1994), or a transient, possibly involving the disintegration of the companion. For a given binary system, the probability of a prompt collision caused by the natal kick (whose direction is isotropically distributed in the frame of the collapsing star) can be estimated as the solid angle subtended by the companion at the position of the collapsing star. For our fiducial population, excluding systems that do not receive a natal kick because of the fallback rescaling, we obtain that roughly 1 out of 10 000 CC events in a binary would result in a prompt collision (i.e., a probability per binary system of about 10⁻⁴). For this estimate, we exclude from the normalization systems merging before the first CC.

8. Astrophysical Implications

8.1. Other possible contributions to the inferred velocity

Throughout this study, we have used velocities calculated in the rest frame of the initial binaries. In reality, the progenitor massive binaries are likely to form in high-mass star forming regions with a certain velocity dispersion, typically smaller than 10 km s⁻¹ (de Bruijne 1999; Steenbrugge et al. 2003; Kiminki & Smith 2018). Moreover, the high-mass star forming regions have been observed to have a systematic lag velocity of about ∼5 km s⁻¹ compared to the rotation of the Galactic disk (Reid et al. 2014).

In an observed sample using the putative parent association to define a frame of reference, the “thermal” velocity of the progenitor binary within the star-forming region will add to the velocity resulting from the disruption. To illustrate how this affects the velocity distribution of Fig. 5, we numerically convolve the distribution for ejected stars more massive than 15 M_⊙ (red) with a Gaussian distribution with full-width-half-maximum of 10 km s⁻¹, assumed to be an upperlimit of the typical velocity dispersion of OB associations (de Bruijne 1999; Steenbrugge et al. 2003; Kiminki & Smith 2018). We note that both the orientation of the thermal velocity of the parent binary, and the velocity acquired by the ejected star are randomly oriented, and the latter depends on the binary inclination and phase at the first CC event.

Figure 12 shows the re-sampled distribution and the result of the convolution with a Gaussian velocity dispersion distribution with dispersion of 10 km s⁻¹, which we consider an over-estimate for OB-associations. While adding a dispersion velocity smears out the peak of the distribution and increases the contribution of the runaways faster than 30 km s⁻¹, the effect is insufficient to reconcile the observed runaway fraction with our results (see Sect. 8.2). The ratio of walkaways to runaways more massive than 15 M_⊙ remains even accounting for the velocity dispersion of the parent association (cf. ℛ₁₅ ∼ 27 before accounting for it).

Fig. 12. The yellow histogram shows the convolution of the distribution of ejected MS stars more massive than 15 M⊙ (red histogram, cf. Fig. 5) with a Gaussian with full-width-half-maximum of 10 km s−1 representing (an upper limit to) the velocity dispersion of a star forming region (shown in the inset plot). In this figure, the red histogram is normalized so that its surface area is one.

For observed samples using the Galactic disk rotation to define a reference frame, the peculiar lag velocity of the high-mass star-forming region could potentially contribute to overestimating the ejection velocity of the stars. If the lag velocity of the parent association is not removed from the measures, it would introduce a systematic shift in one component of the velocity of the star (corresponding to the direction of the peculiar motion of the parent association). This can be modeled with a convolution of our distributions with a Dirac’s δ-distribution with a velocity v_lag = (v_lag, 0, 0) for an appropriate choice of the orientation of the (Cartesian) frame. The direction of the lag velocity is in this case fixed, and only the ejection velocity is randomly oriented.

The convolution with a δ-distribution would also shift the peak of Fig. 5 to higher velocities, but would smear it out much less. Adopting a lag velocity of v_lag ≃ 5 km s⁻¹ (Reid et al. 2014) would barely change the fraction of apparent O-type runaways we predict. To infer an O-type runaway fraction of ∼10%, the true lag velocity of high-mass star forming regions would need to be larger than ≳15 km s⁻¹, and would need not to have been corrected when making inferences from observations.

We note that for O-type runaways with lifetimes of ≲15 Myr, it is likely possible to locate the parent association and correct for its possible peculiar motion. The possible lag velocity of high-mass star forming regions only affects stellar velocities measured using a frame co-rotating with the Galactic disk.

8.2. On the runaway fraction

Our results predict a runaway fraction among O type stars of at best a few percent, in agreement with Eldridge et al. (2011), but in tension with the observational result that ∼10 − 20% of O type stars are runaways (Blaauw 1961; Gies 1987; Stone 1991; Tetzlaff et al. 2011; Maíz Apellániz et al. 2018; Sana et al., in prep.), and the claim that roughly two-thirds of runaways come from the disruption of binaries (Hoogerwerf et al. 2001). The latter has also been challenged by the more recent observations from Jilinski et al. (2010). As a consistency check we can estimate the runaway fraction with:

(6)

where we can assume a binary fraction f_bin = 1, reasonable for the population of O type stars considering initial periods as large as 10^5.5 days, the fraction of binaries that do not merge is 1 − f_mergers ≃ 0.8 (Sana et al. 2012), the disruption fraction is 𝒟 ≃ 0.8 (this study, but also, e.g., Eldridge et al. 2011), the fraction of material accreted, for which we assume f_acc ≃ 0.2 (Packet 1981), and the fraction of ejected companions fast enough to be a runaway which we can estimate¹⁰ from 𝒟 and ℛ₁₅, and is f_{v >  30} ≃ 𝒟/(1 + ℛ₁₅)≃0.04 for our fiducial run. Assuming an IMF slope α = −2.3 we obtain a runaway fraction of f^RW ≃ 0.01, in reasonable agreement with our results.

The tension between our results and the observed runaway fraction suggests that either (i) binaries tend to evolve towards shorter pre-CC orbits, corresponding to higher ejection velocities, but are still disrupted easily by the CC event, (ii) the contribution of dynamical ejections to the runaway population is presently underestimated, (iii) other ejection mechanisms exists and are presently overlooked, (iv) the fast runaways might be easier to detect since they likely move away from their gas-rich birth environment (leading to typically lower extinction, Maíz Apellániz et al. 2018), (v) the velocities of observed O type runaways are overestimated (cf. Sect. 8.1), or any combination of the previous.

Since , with J_orb orbital angular momentum, to shift the peak of the velocity distribution in Fig. 5 from ∼6 km s⁻¹ to v_dis ≳ 30 km s⁻¹, the binaries would need to lose about five times more orbital angular momentum during the evolution before the collapse of the primary. This discrepancy deserves to be revisited once the runaway fraction can be evaluated from the homogeneous Gaia datasets. If the contribution from dynamical ejections can robustly be quantified, the population of runaways might provide new constraints on the mass and angular momentum losses during mass transfer in binaries.

A mechanism to remove angular momentum and mass from the binary presently not included in our population are eruptive mass loss events for which there is growing evidence from early observations of SNe (e.g., Khazov et al. 2016). Such mass loss events might be much more common than previously thought. One possible physical cause for these event are gravity waves excited by late shell burning depositing energy at the base of the envelope (Quataert & Shiode 2012; Fuller 2017; Fuller & Ro 2018). It is possible that the mechanism driving these mass loss events might also result in pre-CC binary interactions which might significantly change the pre-CC orbit. SN-impostor events happening earlier in the evolution might also lead to similar effects on the orbital evolution, which we presently do not model.

Another speculative possibility to increase the angular momentum losses could be magnetic braking. If non conservative mass transfer can simultaneously generate a magnetic field and eject some mass from the system, this might result in torques on the binary much larger than we consider here.

8.3. How to identify walkaway stars

The velocity dispersion of OB stars in our Galaxy is ≲10 km s⁻¹ (e.g., Blaauw 1956; Gies 1987; TT98; Hoogerwerf et al. 2001), and it is typically lower for OB associations (de Bruijne 1999; Steenbrugge et al. 2003; Kiminki & Smith 2018). We emphasize also that the velocity dispersion increases going from samples of OB stars in clusters, in associations, and in the field (Gies 1987). The average velocity of ejected secondaries is typically higher than this value (cf. Table 1), and if considering only the post-binary interaction secondaries (i.e., removing all the ejected secondaries from very wide and non-interacting binaries, cf. dashed distributions in Fig. 5) it increases further, especially for M_dis ≥ 7.5 M_⊙. Therefore, the companions ejected after a mass transfer episode might form a distinct population in the Gaia data. Moreover, regardless of their previous history, the walkaway and runaway stars from binary disruptions do not need to follow the local Galactic rotation curve, and their motion can bring them to locations in which they are more easily observed (e.g., Boubert & Evans 2018; Maíz Apellániz et al. 2018).

Tetzlaff et al. (2011), following Stone (1991), proposed to fit the velocity distribution of young (≲50 Myr) stars within 3 Kpc of the Solar neighborhood using two Maxwellians: one for the “low-velocity” (v ≲ 30 km s⁻¹) and one for the “high-velocity” group. Our findings suggest that ∼10% of the O-stars in the low velocity group might still be the result of binary disruptions.

Because of the binary interactions taking place before the first CC, these stars carry observational signatures which might make them recognizable even if they do not stand out from a kinematic point of view. For instance, if there is a mass transfer phase during the previous binary evolution, the ejected star will be spun up (e.g., Packet 1981; Pols et al. 1991; Boubert & Evans 2018) and possibly chemically polluted with He-rich and/or N-rich material (e.g., Blaauw 1993). Mass transfer might also be responsible for the presence of strong magnetic fields, (e.g., Schneider et al. 2016). These features should in principle leave observable signatures in the spectra of these stars (e.g., Maíz Apellániz et al. 2018). However, the effects of the metallicity gradient in the radial direction of the Galaxy complicates the spectral identification of binary products, as pointed out by McEvoy et al. (2017).

Characterizing the population of stars ejected from isolated binary disruptions is necessary to reduce the number of false positive detections of hypervelocity stars in current and upcoming astrometric catalogues (e.g., Marchetti et al. 2017).

Stars ejected by a successful SN explosion recently (i.e., less time ago than the lifetime of the SN remnant) can be connected to the SN location, and the associated NS can sometimes also be found. Assuming a typical visible lifetime of the SN remnant of τ_SNR ≃ 10⁴ years, the typical distance traveled by an ejected companion is of order D ≡ v_dis × τ_SNR ≃ 0.1 pc. Therefore runaways, but especially walkaways are expected to still reside within the SN remnant while it is still observable. This strategy has been successfully applied to search for runaways by van den Bergh (1980), Guseinov et al. (2005), Tetzlaff et al. (2013, 2014), Dinçel et al. (2015), Boubert et al. (2017b). Kerzendorf et al. (2019) also used the predicted kinematics of disrupted binaries to search for companions of the star that exploded producing the supernova remnant Cas A. Tetzlaff et al. (2011) also suggested that, with precise astrometry, it might be possible to relate the ejected companion to the remnant of the companion up to ∼5 Myr after the SN explosion, provided that the compact remnant of the former companion is visible (e.g., as a pulsar).

The results of Banerjee et al. (2012) and Perets & Šubr (2012) suggest there might be a spatial distinction between dynamically ejected stars and post-binary runaways and walkaways. Dynamical interactions in a cluster are most efficient early in the cluster evolution (e.g., Oh & Kroupa 2016), well before the first stellar CC event happens. Moreover, the median velocity of cluster ejection is higher than the typical walkaway velocity: Banerjee et al. (2012) derived a median velocity of dynamically ejected stars of 50 km s⁻¹ after 3 Myr of evolution. This value is significantly higher than what we find from binary disruptions. The combination of these two effects would suggest that walkaways and runaways from binary disruptions would generally be closer to their parent cluster than dynamically ejected stars. However, this simplistic prediction is complicated by the fact that binaries can eject stars from the parent cluster outskirts, while most dynamical interactions produce runaways from the center of the cluster. Another complication is the possible ejection of stars as tidal tails in a cluster merger process, (e.g., Lucas et al. 2018).

Evans & Massey (2015) and Neugent et al. (2018) have reported the first discovery of extraGalactic post-MS massive runaways. They reported the observation of a red super giant in M31 and a yellow super giant in the SMC, respectively. They inferred peculiar space velocities larger than ∼400 km s⁻¹ and ∼150 km s⁻¹ for these stars, and proposed a binary origin at least for the yellow supergiant in the SMC. Our fiducial population can hardly reach such high peculiar velocities, unless the velocity of the system as a whole binary itself was already high. Even if these stars were ejected from a massive binary during their main sequence, velocities in excess of a few hundred km s⁻¹ are only marginally reached in our simulations (cf. Figs. 5, 8 and C.2). If our present understanding or binary physics is correct, this discrepancy might suggest a different ejection mechanism (or combination of mechanisms) for these two stars. Evans & Massey (2015) also report that vast majority of O/B stars and red super giant stars in the Milky-Way have radial velocities lower than ∼70 km s⁻¹, in very good agreement with our results.

8.4. X-ray binaries and gravitational wave sources

Binaries remaining bound after the CC of the primary are a minority. Selecting from our results only bound systems consisting of a BH and a MS star, we obtain the systemic velocity distributions for BH binaries with different BH kick scenarios (see right panel of Fig. 11). Gaia will give astrometric constraints on the 19 known BH X-ray binary in the Galaxy: if the effect of the Galactic potential can be singled out from these data, the comparison with our distribution might shed light onto the typical amplitude of BH kicks (e.g., Fragos et al. 2009, 2013; Repetto et al. 2012; Mandel 2016).

The binaries surviving the first CC can become X-ray sources if their separation is short enough for the compact object to accrete during the subsequent evolution. For all our natal kick assumptions, the post-CC separation peaks at ∼300 R_⊙, and most bound system have separations smaller than ∼500 R_⊙, indicating that the majority of the systems remaining bound will become X-ray sources, and potentially even go through a common envelope evolution phase because of the large mass ratio between the compact object and the secondary star (which has accreted mass, cf. Sect. 3).

However, observed X-ray sources more likely have relatively short orbital periods, or, in other words, the energy released at the CC of the first star did not widen significantly the binary and instead contributed mostly to the kinetic energy of the system. This could mean that the observed sample of Galactic BH X-ray binaries is biased towards the high velocity tail of our distributions.

System that might become X-ray sources have systemic velocities of a few tens of km s⁻¹, and since the secondary star is also an accretor and rejuvenated, they are able to travel a distance comparable to ⟨L_walk⟩. This can also change their contribution as feedback engines (e.g., Justham & Schawinski 2012; Fragos et al. 2013).

If a bound system is not disrupted by the first CC and does not merge during the evolution of the secondary star, it is less likely to be disrupted by the second CC (e.g., Stevenson et al. 2017; O’Shaughnessy et al. 2017), and therefore more likely to become a gravitational wave source (if the separation is small enough). This is because the impact of the kick on the orbit scales as v_k/v_orb, and the orbit of bound systems is likely to shrink during the evolution after the first CC: the donor star needs to become less massive than the compact object for the system to widen as described in Sects. 3 and 4. We also emphasize that typically RLOF has a larger impact than stellar winds on the orbital separation.

8.5. Runaways and walkaways as feedback engines

The large majority of binary systems are disrupted at the first CC (with a few notable exceptions, see also Sect. 6): systems that can become X-ray sources and/or gravitational wave sources are the exception rather than the rule. The constraints on binary evolution coming from the much more common evolutionary path should be considered together with those coming from the more rare channels.

Because of their motion, both walkaways and runaways can reach significant distances from their birth location, see Table 1 and Fig. 6. The typical distances traveled by the ejected secondaries (cf. Fig. 6) are sufficient to get massive stars out of their parent cluster or association, but only the fastest runaways would get beyond the thick disk in our Galaxy, especially if considering the distribution of relative inclinations between the trajectories of the stars and the disk itself. As noted by Boubert & Evans (2018), this introduces a bias in observed samples that favors the detection of the fastest runaways.

The motion of the ejected stars effectively spreads the massive stars and their (possible) SN explosions in a larger volume, enhancing their impact on the composition of galaxies (more efficient mixing of nuclear yields), on the ionization of their surroundings (enhanced escape fraction, e.g. Conroy & Kratter 2012; Kimm & Cen 2014; Ma et al. 2016), but it might possibly decrease the efficiency of their mechanical feedback because of the lower density of the surrounding gas. The large number of walkaways formed per each runaway ℛ suggests that walkaway stars might change massive stars feedback more significantly than the faster runaways, although they cannot reach as far out as the latter because of their lower velocity.

As a consequence of the typical evolutionary scenario described in Sect. 3, the first CC event in a typical binary system will result either in a stripped SN or a failed SN with little or no electromagnetic signature. The second CC will instead be from the accretor star, which could again result in a SNIb/Ic (if the star becomes a Wolf-Rayet because of its own stellar wind), but, because of the IMF and the possibly non-fully conservative accretion rate, the secondary is more likely to remain H-rich until its CC. As already pointed out by Eldridge et al. (2011), this suggests that stripped SNe should preferentially happen in star forming regions, while type II SN should show a larger spatial spread (because of the motion of the ejected star).

In this study, we did not consider the effects of the background cluster and/or Galactic potential on the motion of the ejected secondaries. These can in principle be used as test particles to probe the potential itself, and consequently allow for testing the local and/or Galactic dark matter distribution (e.g., Rossi et al. 2014; Marchetti et al. 2017).

While runaways have peculiar kinematic characteristics that make them recognizable, walkaway stars might more easily be mistaken for genuine single stars. Especially for masses larger than 7.5 M_⊙ and velocities larger than v_dis ≳ 20 km s⁻¹, the walkaways and runaways accrete mass before being ejected. Therefore, the disruption of binaries can pollute observed samples of present-day single stars with binary evolution by-products populating the field of the Galaxy (Gvaramadze et al. 2012; de Mink et al. 2014).

As an example, the runaway star ζ Puppis has been considered to be a canonical O-type star (e.g., Ramiaramanantsoa et al. 2018), and it was used to calibrate the free parameters in the Castor et al. (1975) (CAK) theory of stellar winds (Pauldrach et al. 1994) widely applied to model single stars. However, this object is likely to be an accretor ejected from a binary, (van Rensbergen et al. 1996; Ramiaramanantsoa et al. 2018).

Massive walkaway stars ejected from a binary might also appear in isolation. The origin of massive stars observed in isolation is an open problem: did they form in isolation? or could they possibly have reached that location, if ejected from a binary? Bestenlehner et al. (2011) identify an isolated ∼150 M_⊙ star (VFTS682) in the 30 Doradus field. Such a massive object is unlikely to be the result of a binary disruption, but might have been ejected from the massive young cluster R136. Bestenlehner et al. (2011) suggest that if this star was not ejected from this cluster, it might be a direct indication of isolated star formation. Walkaway massive stars can pose similar problems if not recognized as rejuvenated binary products.

This might also connect to the evolutionary path of LBV stars. These are massive stars that experience extreme and yet unexplained outbursts of mass loss. In the classical picture of massive stellar evolution (so-called Conti scenario), LBVs are an intermediate evolutionary stage between O-type and Wolf-Rayet stars, during which the bulk of the H-rich envelope is lost (e.g., Conti 1975; Maeder & Conti 1994; Maeder 1996). However, Smith & Tombleson (2015) pointed out that LBV and LBV-candidates are farther from the nearest O-type star than Wolf-Rayet stars, and suggested that LBVs might be rejuvenated products of binary evolution¹¹ (accretors or merger products). The follow-up study of Aghakhanloo et al. (2017) suggested that partial rejuvenation during mass transfer coupled to a walkaway (or runaway) velocity might be sufficient to explain the isolation observed by Smith & Tombleson (2015).

9. Summary and conclusion

We have carried out a suite of numerical simulations to provide predictions for the kinematics and properties of massive stars ejected from binary systems disrupted by the supernova explosion of the companion star. Our aim with this study was three-fold: (i) identify which theoretical predictions are robust against model uncertainties and discuss their astrophysical implications (ii) investigate theoretically which uncertain physical processes affect the kinematics of the population of unbound stars in such a way that their imprints would be observable and (iii) provide a framework of models that can be used for comparison against observations or as input for other simulations.

We summarize our main findings below. When we quote percentages or fractions, the error bars given indicate the maximum variations that we encountered in our model variations for different physical assumptions (excluding the unphysical variation with zero natal kicks for all newly born compact objects).

• Nearly a quarter, %, of all binary systems with at least one star more massive initially than 7.5 M_⊙ merges prior to the explosion of the first star (see Fig. 4 and Table 1).

• The large majority of binary systems are disrupted at the moment of core collapse of the first star, in our simulations, consistent with earlier studies. This prediction is robust against variations in the treatment of the evolutionary and interaction processes, but it is sensitive to the choices concerning the natal kicks of neutron stars and black holes.

• The disruption of the binary system produces an unbound stellar companion. Remarkably, we find that these events only rarely produce runaway stars (i.e., stars with peculiar motions larger than 30 km s⁻¹). The velocity distribution peaks at ∼6 km s⁻¹, with a 90th percentile around 20 km s⁻¹ well below the typical threshold for runaway stars.

• The slow-moving unbound companions, which we call walkaway stars, outnumber the runaway stars by at least an order of magnitude in nearly all our models. This is a robust outcome of the models and is primarily due to the reversal of the mass ratio prior to the explosion as we show numerically and analytically. This finding has been noted in earlier studies (e.g., De Donder et al. 1997; Eldridge et al. 2011), but we believe this result (and its potential implications) is underappreciated.

• We also produce runaway stars but we rarely find them to be faster than 60 km s⁻¹, which corresponds to 98.9th percentile of our fiducial population. This is much lower than the theoretical maximum set by the finite size of the stars and their orbit, which is in the range 400–1000 km s⁻¹ depending on their masses and evolutionary state. We expect that stars with larger peculiar velocities result from dynamical ejections from clusters, capture from another Galaxy, or (most exotically) interaction with an intermediate mass or supermassive black hole.

• Isolated binary evolution produces a runaway fraction among stars more massive than 15 M_⊙ between 0.2 and 2.6%. None of the combination of parameters that we consider produces a runaway fraction approaching the observed value of ∼10% for O-type stars. This result is consistent with previously published theoretical studies (see also Appendix A), and might indicate that (i) the observed fraction of O-type stars that are runaways is overestimated, possibly because of biases favoring the observation of fast and isolated massive stars or because of the difficulties in defining the frame of reference in which to measure the velocity, or (ii) contrary to previous claims, it is not correct that the majority of observed O-type runaways come from the disruption of binaries, or (iii) massive binary evolution models are lacking some physical process allowing them to produce short pre-collapse orbits which would result in faster ejection velocities if disrupted.

• We find a mild trend for runaway stars to be faster and more common at lower metallicity, caused by the fact that stars are on average smaller and the reduced effects of stellar winds. Gaia should in principle be able to test these predictions by identifying the fastest runaway stars in the Magellanic clouds. In our most metal poor simulations for Z = 0.0002, which may have relevance for the earliest stellar populations, we find five times more runaway stars with masses larger than 7.5 M_⊙, if we keep all other assumptions for the initial distributions the same. The sensitivity of our results to the assumed metallicity also suggests that binary evolution might create systematic trends in the stellar feedback with metallicity and galaxy size.

• Both runaways and walkaways typically accrete mass from their companion prior to the disruption, This is especially true for those more massive than 7.5 M_⊙ at the time they are ejected. They typically also gain angular momentum and nuclearly processed material in the process. These features might make walkaway stars recognizable as binary products even if they do not stand out from a kinematic point of view. Their further evolution and structure may also differ from that of non-rotating single stars, as has been speculated in the context of the LBV phenomenon.

• We suggest that the high mass tail of the runaway mass function could provide insight on the formation mechanism of black holes from the core-collapse of the primary star, and in particular on whether these receive a “velocity” or a “momentum” kick. So far, studies have focused on the surviving X-ray binaries, but we propose using their unbound counterparts since they are more common and should be identified by Gaia.

• By ejecting secondaries, the disruption of binaries can enhance the role of massive stars as feedback engines. The average distance they travel, calculated neglecting any external potential, is of order ≳100 pc, and up to ≳500 pc if considering only the faster runaways. These values suggest that both runaways and walkaways can exit the dust cloud embedding their birth location, thus changing the escape fraction for their ionizing radiation. If they are massive enough, the ejected star will explode far from other stars and gas overdensities, with potential implications for the driving of turbulence and star formation.

The results for each parameter variation are available online¹² and at CDS link.

Acknowledgments

For the many helpful discussions we acknowledge F. Broekgaarden, A. de Koter, D. Hendriks, L. Kaper, N. Langer, I. Mandel, C. J. Neijssel, Ph. Podsiadlowski, K. Sakar, F. N. Schneider, S. N. Shore, E. van den Heuvel, and A. Vigna-Gomez. MR is grateful to the KITP for the visiting graduate fellowship during which a large portion of this study was carried out. This research was supported in part by the National Science Foundation under Grant No. NSF PHY11-25915. SdM has received funding under the European Unions Horizon 2020 research and innovation programme from the European Research Council (ERC) (Grant agreement No. 715063). RGI thanks the Science and Technology Facilities Council (STFC) for funding his Rutherford fellowship under grant ST/L003910/1.

References

Appendix A: Comparison to selected previous population studies

Various earlier studies have simulated the production of binaries containing a compact object and/or ejecting runaway stars. We briefly discuss a selection that focuses on estimates for unbound companions and comment briefly on how our findings compare.

De Donder et al. (1997) presented an extensive binary evolution study discussing predictions for the O-type runaways as well as the systems that remain bound. They consider a variety of assumptions for the initial distributions and the uncertain physical parameters concerning the efficiency of mass transfer, angular momentum loss and the common envelope ejection. Generally, we find that our results agree well. They find that 16–23% of systems remain bound after the SN explosion of the primary, which is consistent with our findings . They find that between 6% and 27% of the O-type are unbound former companions that are now single, which is also agreement with our walkaway fraction (one but last column in their Table 1).

The find that 2–7% of the O-type stars have velocities large than 30 km s⁻¹ (final column in their Table 1), which is slightly larger but still consistent with what we find. We expect that this difference may be in part the result of a difference the treatment of mass transfer and angular momentum loss, although the combined effect of further differences in our assumptions will also contribute. Their default assumption is that a fixed fraction β_RLOF = 0.5 of the material transferred during Roche-lobe overflow leaves the system through the outer L2 Lagrangian point forming a ring around the binary system, which leads to larger angular momentum loss, shrinking the orbit further. This is somewhat similar to our model variation where we assumed γ_RLOF = γ_disk although our results cannot be compared one to one, since we adopt a different, physically motivate, assumption for the mass transfer efficiency. In our simulations, we find that this assumption is not increasing the number of runaway stars, because the fraction of system that merges increases too.

Eldridge et al. (2011) perform population synthesis simulations with a detailed stellar evolutionary evolutionary code to investigate O and early B type runaway stars. Their primary aim is the predict the spatial distribution of different types of CCSN and gamma-ray bursts. They find a disruption fraction, 𝒟 = 80%, in good agreement with our results, . In their study, they use the term runaway to refer to all unbound companions with velocities larger than 5 km s⁻¹, which encloses the large majority of what we refer to as walkaway stars, but they also quote estimates for stars faster than 30 km s⁻¹ which corresponds to our definition of runaway stars. They predict a walkaway (runaway) fraction of 2.2% (0.5%) for O type stars in their simulations for Z = 0.02. For comparison, for stars more massive than 15 M_⊙, which roughly corresponds to O type stars, we estimate a walkaway (runaway) fraction of % (%). Our runaway fraction agrees very well. We find a somewhat larger walkaway fraction, but we consider this a fairly good agreement, given the uncertainties and differences in definitions that we have adopted. Also the distance they estimate that ejected stars can reach is in good agreement with our predictions and they also note that the majority unbound companions accrete mass from their companions before the first core-collapse. Our predictions for the velocities of bound systems are slightly lower compared to Eldridge et al. (2011): the vast majority of our bound post-core-collapse systems are slower than 30 km s⁻¹, and there is a high-velocity tail barely extending beyond about 100 km s⁻¹.

Very recently, Boubert & Evans (2018) published a study investigating the hypothesis that Be stars are products of mass transfer in binary systems (e.g., Pols et al. 1991; de Mink et al. 2013). They compare the kinematics of a flux-limited sample of Galactic Be stars with binary population synthesis simulations. These simulations were obtained with a different version of the binary evolutionary code that we use. Generally, our findings are in agreement, despite the minor differences in the model assumptions. They also find a large fraction of unbound companions, many of which are rapidly rotating and moving at velocities slower than 30 km s⁻¹ (e.g., their Figs. 6 and 7b), in good agreement with our Fig. 7. They also find that the natal kick distribution does not greatly affect the resulting velocity distribution and/or the runaway fraction and that the maximum distance traveled by ejected Be stars is likely to be smaller than the vertical scale height of the thin disk, in agreement with our findings.

We find some disagreement in the provided explanations of the theoretical results. For example, the authors state in their Sect. 3.2 that that mass transfer shrinks the orbit and accelerates the secondary star. While this the case initially, upon the onset of mass transfer, we find that the orbit generally widens after the reversal of the mass ratio. The widening and the inversion of the mass ratio both slow down the orbital velocity as we verified both analytically Sect. 4 and by detailed inspection of representative example systems, e.g., Sect. 3. In their Sect. 3.3, they state that “whether a binary is disrupted by a supernova is principally determined by whether the primary loses more than half its mass (Blaauw 1961), and the kick on the compact object is only a second order effect”. We find instead that mass loss during the explosion is rarely responsible for unbinding of a binary system that can produce a runaway star. The amount of mass lost needs to exceed half of the total mass of the system (and not half of the primary star), which is rarely achieved in our simulations. This is because the CC progenitor loses most of its mass during the preceding mass transfer phase. We find that the Blaauw kick due to rapid mass loss is only important of initially very wide binaries, in which the two stars have not exchanged mass prior to CC. This difference may be in part due to the differences in the assumptions for the range of initial orbital periods. The authors consider systems with initial orbital periods up to 10¹⁰ days in their simulations, which means that the majority of their progenitors should effectively evolve as single stars. We consider instead systems up to 10^5.5 days, which is more appropriate for the more massive progenitors that we are interested in.

Many further studies investigated the populations of binaries in the context of the formation of X-ray binaries and binary neutron stars and black holes. Providing a complete overview and detailed comparison is beyond our present scope. We discuss below a comparison to a limited set of studies.

Brandt & Podsiadlowski (1995) focus on binary systems remaining bound at the first CC, with the aim of understanding the effects of natal kicks on X-ray binaries. They estimate a disruption fraction in their calculations of 𝒟 ≃ 73 − 81% (depending on the companion mass), in good agreement with our results.

Kalogera (1996) presented a similar analytic study of the effects of natal kicks on the systemic velocities of X-ray binaries. Our population synthesis results agree in predicting systemic velocities of bound post-CC systems generally lower than the pre-explosion orbital velocity, except with (corresponding to large ξ in the notation used by Kalogera 1996).

Fryer et al. (1998) investigated the impact of NS natal kicks on the formation of NS X-ray binaries and NS-NS binaries, but also present results for the unbound companions. Assuming a bimodal kick distribution, they find that most O/B-type ejected companion move slower than about 50 km s⁻¹, consistent with our findings.

Dray et al. (2005) focus on high-mass runaways, which become Wolf-Rayet stars during or before their post-disruption evolution. They argue in favor of significant BH kick amplitudes to explain the rarity of BH-WR binaries and the observed velocity distribution of WR runaways. We reach similar conclusions based on our simulations for different assumptions for the BH kicks. The systemic velocities we find for bound systems are also in reasonable agreement with those found by Dray et al. (2005)

Recently, Tauris et al. (2017) presented a detailed study of the evolutionary processes leading to the formation of NS NS binary systems. They also find that the majority of binaries hosting a NS after the first CC have systemic velocities smaller than ∼30 km s⁻¹, in good agreement with the results shown in the left panel of Fig. 11.

Appendix B: Output files

The outcome of our population synthesis calculations are available online¹³ and at CDS link. Each file corresponds to one parameter variation (see Sects. 2 and 6), and logs the following information for each binary system where a star goes CC¹⁴.

• primary ZAMS, pre-CC, and post-CC (corresponding to the NS or BH mass if the primary is the star collapsing) masses in M_⊙ units: M1zams, M1preCC, M1postCC;

• secondary ZAMS, pre-CC, and post-CC masses in M_⊙ units: M2zams, M2preCC, M2postCC;

• fallback fraction f_b for each star (set to zero for the star that is not collapsing): fb1, fb2;

• evolutionary stage before and after CC (stellar types listed according to Hurley et al. 2000): type1preCC, type1postCC, type2preCC, type2postCC;

• post-CC velocities of the 2 stars in the original frame, in km s⁻¹: v1postCC, v2postCC;

• pre-CC and post-CC eccentricity (the latter is -1 for mergers and disrupted systems): e_preCC, e_postCC;

• pre-CC and post-CC separation in R_⊙ units (the latter is set to 0 for mergers and disrupted systems): a_preCC, a_postCC;

• ZAMS, pre-CC, and post-CC periods in days: Pzams, PpreCC, PpostCC;

• kick amplitude in km s⁻¹: v_kick;

• kick direction, with θ angle between the collapsing star orbital velocity and the kick (see also TT98 for notation): theta, phi;

• systemic velocity in km s⁻¹: v_sys;

• age of the system at the time of CC in Myr: t_explosion;

• time left in the current evolutionary stage in Myr: t_remaining;

• time spent by the system with at least one star more massive than 15 M_⊙: duration_*;

• system probability (see below): Prob.

The system probability corresponds to the hyper-volume of the initial parameter space (, q^ZAMS, P^ZAMS) represented by each binary system in our model grid. In other words, the probability of each system is the statistical weight of the system seen as sampling point for the initial distributions. To construct distributions of the output quantities (e.g., those presented in Figs. 5–11), the properties of each system should be weigh with the corresponding probability. Similarly, the mean value of a quantity ⟨x⟩ (e.g., ⟨v⟩ in Table 1) should be calculated using:

(B.1)

where P is the probability, and the index i runs over all the binary systems in a population.

To calculate the average distance traveled by stars ejected by the binary disruption, we use:

(B.2)

where the second term in squared brackets accounts, albeit in a simplified way, for the helium core burning duration of the rejuvenated star.

Appendix C: Observable velocity distribution

Figure C.1 shows the velocity distribution of MS stars ejected by the disruption of binaries that can be directly compared to observations provided that (i) the contribution of dynamical ejection can be separated in the observed sample and (ii) the effects of the Galactic potential can be neglected, i.e. effectively each ejected star moves in a straight line at constant velocity for the remaining duration of its MS. This is the same information of Fig. 5, but each bin is populated considering also the remaining MS lifetime of the ejected star (τ_MS).

Fig. C.1. Velocity distribution of ejected stars, including the finite MS lifetime to populate the bins (see also Fig. 5).

Fig. C.2. Same as Fig. 5, but using a logarithmic scale for the velocity. The use of a logarithmic scale allows for the display of a wider range of velocities. A minor peak in the grey histogram can be seen between 100 ≲ vdis/km s−1 ≲ 400, but is absent in the histograms for massive ejected stars. Such high ejection velocities are reached through a common envelope evolution without accreting mass.

Appendix D: Pre-collapse distributions

We present in this appendix the pre-CC distribution in separation, mass of the collapsing star, and mass of the companions for all the binaries with a MS companion to the collapsing star in our fiducial simulation. Similar distributions can be derived for all our parameter variations from the data files that will be made available. These distributions can inform studies of the interaction of the SN shock with the companion star (e.g., Wheeler et al. 1975; Liu et al. 2015; Rimoldi et al. 2016; Hirai et al. 2018).

We show in Fig. D.1 the pre-CC separation distribution. The colors indicate indicate the minimum mass of the MS companion (i.e., not of the collapsing star). Roughly speaking, the two peaks shown in Fig. D.1 correspond to the orbital widening due to conservative (case A and early case B) and non-conservative (late case B and case C) mass transfer phase. For all pre-CC separations shorter than 10³ R_⊙, the dashed and solid histograms coincide, indicating that all these binaries have experienced a direct interaction previously during the evolution.

Figure D.2 shows the mass distribution for the exploding star and the MS companion (M_CC, and M₂, respectively), at the pre-CC stage. The combination of the distributions shown in Figs. D.2 and D.1, together with the effects of the natal kick distribution results in the ejection velocities in Fig. 5 which is our main result.

Fig. D.1. Pre-CC separation distribution for binaries with a collapsing star and a main sequence companion. Colors indicate the pre-CC mass of the MS companion according to the legend. Dashed histograms indicate post-interaction (RLOF or common envelope) binaries. Top panel: corresponding cumulative distributions.

Fig. D.2. Pre-CC mass distribution for the exploding star (MCC) and the companion (M2). We plot all systems where the companion is a MS star at the time of the explosion, regardless of whether the binary is disrupted or not. Dashed histograms indicate post-interaction (RLOF or common envelope) binaries. Top panel: corresponding cumulative distributions.

All Tables

All Figures

Fig. 1. Schematic depiction of the evolution of a close massive binary through stable mass transfer. The evolution forks at the first core-collapse event in the system: the vast majority of the systems are disrupted and each produces a runaway or walkaway star which can travel far from its birth location. The systems remaining bound are possible progenitors of X-ray binaries and gravitational wave sources. The fraction of binaries disrupted 𝒟 comes from the simulations presented in this study.

In the text

Fig. 2. Evolution of an example binary system with = 20 M⊙, = 15 M⊙ with PZAMS = 100 days, described in Sect. 3 (see also the cartoon in Fig. 1). Top panel: primary, secondary, and total mass. Bottom panel: orbital velocity of the secondary (left y-axis) and period of the binary (right y-axis). The labels indicate the following phases: A. the MS of the primary, B. case B RLOF, emphasized by the vertical dashed line, C. the binary is made of a He star (possibly looking like a Wolf-Rayet star) and a rejuvenated MS secondary; D. (and the vertical dot-dashed line) marks the explosion of the primary in a SNIb/Ic which unbinds the NS and its companion. During the phase E1. the secondary is a slow moving walkaway star, traveling for a duration of about half the main sequence lifetime of the primary.

In the text

Fig. 3. Analytic estimates of the final-to-ZAMS ratio of the specific orbital angular momentum h (left panel), the orbital separation a (center panel), and orbital velocity of the secondary v2 (right panel), assuming qZAMS = 0.75 and μenv = 0.75, which are representative for the example system of Sect. 3. Blue colors correspond to a decrease, red to an increase. The mass transfer efficiency βRLOF and the angular momentum parameter γRLOF are assumed constant throughout the evolution.

In the text

	Fig. 4. Overview of the binary evolution scenarios up to the first CC event. The branching ratios shown are from our fiducial simulation, and we highlight in red the disruption fraction 𝒟. The errors on each fraction exclude the run without SN kicks (σkick = 0 km s−1), which produces an unrealistically low disruption fraction (cf. Table 1 and Sect. 6).
In the text

Fig. 5. Velocity distribution of MS stars ejected from a binary system at the time of the first CC. Top panel: corresponding cumulative distributions. The walkaways are one order of magnitude more numerous than the runaways. The gray, blue, and red histograms show all MS secondaries, only MS stars more massive than 7.5 M⊙ (i.e., roughly those that might experience CC), and only MS stars more massive than 15 M⊙ (i.e., roughly the O-type stars), respectively. The mass considered here is taken right after the CC of the primary star. The dashed lines show the distribution of walkaways and runaways that have gone through RLOF (or common envelope evolution) before being ejected. Almost all massive runaways and walkaways, and the majority of the ejected stars at large velocities, have gone through RLOF. The inset plot magnifies the runaway regime vdis ≥ 30 km s−1. See also Fig. C.2 for a wider velocity range, Fig. C.1 for a figure accounting for the remaining lifetime as walkaway or runaway star, and Fig. 12 for an example of how the velocity dispersion of a star forming region affects these distributions.

In the text

	Fig. 6. Distribution of the maximum distance L that ejected MS secondaries can reach, neglecting any gravitational potential. The blue and red histograms are for ejected stars more massive than 7.5 M⊙ and 15 M⊙, respectively. The vertical dashed lines mark the mean values of the distributions.
In the text

Fig. 7. Equatorial velocity at the time of ejection for massive (Mdis ≥ 7.5 M⊙) companions as a function of the ejection velocity. Brighter colors indicate location of the parameter space more populated by our fiducial simulation. All the massive runaway stars spin with an equatorial velocity of ∼500 km s−1, close to breakup rotation, since they have accreted mass from their companion before being ejected. The spread is due to wind spin down before the binary disruption. We do not include the effect of projecting on the line of sight, or post-ejection wind spin down in this plot. The vertical dot-dashed line marks the threshold to define runaway stars.

In the text

Fig. 8. Velocity distribution of companions ejected after their MS. Helium stars and white dwarfs (WD) can have much shorter pre-CC orbits because of their smaller radii, and thus they can reach much higher ejection velocities. Conversely, the number of non-degenerate, hydrogen-rich post-MS stars is dominated by red supergiants with very large radii, which can only exist at large separation from the exploding companion, and thus typically have much slower velocities. Top panel: corresponding cumulative distributions.

In the text

Fig. 9. Mass function of massive (Mdis ≥ 7.5 M⊙), runaway (faster than ≥30 km s−1) stars depends on the assumptions for the BH kick. Bottom panel: our fiducial simulation with σkick = 265 km s−1, including the fallback downscaling of the kick both for NSs and BHs. The fallback is an important effect only for BHs. Central panel: mass function for intermediate BH kicks (σkick = 100 km s−1, no fallback downscaling) and our fiducial kick for NSs. Top panel: distribution with our fiducial kick amplitude, but without fallback downscaling (i.e., BH velocity kick). The dashed lines indicate the 90th percentiles of the mass distributions.

In the text

Fig. 10. Mass-velocity correlation for the three BH kick variations discussed in Sect. 6.2. The brighter colors indicate where the probability per system is higher on this plane, the colors have a logarithmic scale. Left panel: our fiducial distribution, central panel: case for intermediate BH kicks (σkick = 100 km s−1 and no fallback), and right panel: large BH kicks (same as for NSs). The NS kick is effectively the same in all panels, since the amount of fallback is very small for SNe producing a NS.

In the text

Fig. 11. Normalized distribution of systemic velocities of binaries with a MS star remaining bound after the CC of the first star. The distributions in both panels are normalized to unity. Left panel: distribution for NS companions. Right panel: only binaries with a BH companion. The blue histogram represents our fiducial simulation. The orange solid line gives the distribution when the kick is not rescaled because of fallback, which changes the BH kicks, but not the NS kicks. The red distribution corresponds to a double Maxwellian kick distribution, which changes the NS kicks but not the BH kicks. Top panels: corresponding cumulative distributions.

In the text

	Fig. 12. The yellow histogram shows the convolution of the distribution of ejected MS stars more massive than 15 M⊙ (red histogram, cf. Fig. 5) with a Gaussian with full-width-half-maximum of 10 km s−1 representing (an upper limit to) the velocity dispersion of a star forming region (shown in the inset plot). In this figure, the red histogram is normalized so that its surface area is one.
In the text

	Fig. C.1. Velocity distribution of ejected stars, including the finite MS lifetime to populate the bins (see also Fig. 5).
In the text

Fig. C.2. Same as Fig. 5, but using a logarithmic scale for the velocity. The use of a logarithmic scale allows for the display of a wider range of velocities. A minor peak in the grey histogram can be seen between 100 ≲ vdis/km s−1 ≲ 400, but is absent in the histograms for massive ejected stars. Such high ejection velocities are reached through a common envelope evolution without accreting mass.

In the text

	Fig. D.1. Pre-CC separation distribution for binaries with a collapsing star and a main sequence companion. Colors indicate the pre-CC mass of the MS companion according to the legend. Dashed histograms indicate post-interaction (RLOF or common envelope) binaries. Top panel: corresponding cumulative distributions.
In the text

	Fig. D.2. Pre-CC mass distribution for the exploding star (MCC) and the companion (M2). We plot all systems where the companion is a MS star at the time of the explosion, regardless of whether the binary is disrupted or not. Dashed histograms indicate post-interaction (RLOF or common envelope) binaries. Top panel: corresponding cumulative distributions.
In the text