11institutetext: Astronomical Institute Anton Pannekoek, University of
Amsterdam, 1098 XH Amsterdam, The Netherlands 22institutetext: School of
Astronomy & Space Science, University of the Chinese Academy of Sciences,
Beijing 100012, China 33institutetext: National Astronomical Observatories,
Chinese Academy of Sciences, Beijing 100012, China 44institutetext:
Astrophysics Research Group, Faculty of Engineering and Physical Sciences,
University of Surrey, Guildford, Surrey, GU27XH, United Kingdom
55institutetext: Institute of Astronomy, University of Cambridge, Madingley
Road, Cambridge CB30HA, UK 66institutetext: Institute of Astronomy, KU Leuven,
Celestijnenlaan 200 D, B-3001 Leuven, Belgium

# Massive Runaway and Walkaway Stars

A study of the kinematical imprints of the physical processes
governing the evolution and explosion of their binary progenitors

M. Renzo 11 E. Zapartas 11 S. E. de Mink 11 Y. Götberg 11 S. Justham 2233 R.
J. Farmer 11 R. G. Izzard 4455 S. Toonen 11 H. Sana 66

We perform an extensive numerical study of the evolution of massive binary
systems to predict the peculiar velocities that stars obtain when their
companion collapses and disrupts the system. Our aim is to (i) identify which
predictions are robust against model uncertainties and assess their
implications, (ii) investigate which physical processes leave a clear imprint
and may therefore be constrained observationally and (iii) provide a suite of
publicly available model predictions, to allow for the use of kinematic
constraints from the _Gaia_ mission. We find that $22_{-8}^{+26}$% of all
massive binary systems merge prior to the first core-collapse in the system.
Of the remainder, $86_{-9}^{+11}\%$ become unbound because of the core-
collapse. Remarkably, this rarely produce runaway stars (observationally
defined as stars with velocities above $30\,{\mathrm{km\ s^{-1}}}$). These are
outnumbered by more than an order of magnitude by slower unbound companions,
or “walkaway stars”. This is a robust outcome of our simulations and is due to
the reversal of the mass ratio prior to the explosion and widening of the
orbit, as we show analytically and numerically. For stars more massive than
$15\,M_{\odot}$, we estimate that $10^{+5}_{-8}$% are walkaways and only
$0.5^{+1.0}_{-0.4}$% are runaways, nearly all of which have accreted mass from
their companion. Our findings are consistent with earlier studies, however the
low runaway fraction we find is in tension with observed fractions of about
10%. Thus, astrometric data on presently single massive stars can potentially
constrain the physics of massive binary evolution. Finally, we show that the
high end of the mass distributions of runaway stars is very sensitive to the
assumed black hole natal kicks and propose this as a potentially stringent
test for the explosion mechanism. We also discuss companions remaining bound
which can evolve into X-ray and gravitational wave sources.

###### Key Words.:

stars: kinematics, binaries, massive, supernovae

††offprints: M. Renzo<EMAIL_ADDRESS>
##  1 Introduction

Stars with initial mass larger than about $7.5\,M_{\odot}$ are the progenitors
of black holes (BH) and neutron stars (NS). These stars play an important role
in shaping galaxies through their radiative, chemical, 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 is 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\gtrsim 30\,\mathrm{km\ s^{-1}}$ (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\,\mathrm{km\ s^{-1}}$ 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 is 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). Fig. 1 sketches the typical evolution of a massive binary
system. Most binaries are disrupted at the time of the first core-collapse.

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 ($\gtrsim$$10^{3}\,\mathrm{km\ s^{-1}}$, 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.

Figure 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 is 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
$\mathcal{D}$ comes from the simulations presented in this study.

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., 2016b). 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 et al., 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 Sec. 2).

We first provide insight in our simulations by presenting the results for an
individual system in Sec. 3. We present in Sec. 4 an analytic estimate of the
typical velocity of ejected companions expected from the disruption of
binaries. Then, in Sec. 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 Sec. 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. (2017b); 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 ejects a
slow-moving walkaway star ($v<30\,\mathrm{km\ s^{-1}}$), rather than a faster
runaway. This result is robust against variations of uncertain parameters in
the model, as we discuss in Sec. 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 Sec. 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 will provide our numerical
results upon publication.

##  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; de
Mink et al., 2013; Schneider et al., 2015; Izzard et al., 2018). 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 Sec. 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 primary111Throughout this study, we define the
primary star to be the initially more massive star, even if it becomes the
less massive star in the binary during the evolution. $M_{1}^{\mathrm{ZAMS}}$,
initial mass ratio $q\stackrel{{\scriptstyle\mathrm{def}}}{{=}}M_{2}/M_{1}$,
and initial orbital period $P^{\mathrm{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_{\odot}$ for our fiducial assumptions, assuming there are no
binaries for $M_{1}^{\mathrm{ZAMS}}\leq 2\,M_{\odot}$, see below).

We select $N_{M_{1}}$ primary stars with ZAMS mass $M_{1}^{\mathrm{ZAMS}}$ at
logarithmically spaced intervals in the range $7.5\,M_{\odot}\leq
M_{1}^{\mathrm{ZAMS}}\leq 100\,M_{\odot}$. We weigh each primary star with an
initial mass function (IMF) with slope $\alpha=-2.3$ (Kroupa, 2001). In our
model variations, we explore values of $\alpha$ of -1.9 (Schneider et al.,
2018) and -3.

For each primary star of mass $M_{1}^{\mathrm{ZAMS}}$, we select $N_{q}$
secondaries with mass $M_{2}^{\mathrm{ZAMS}}=qM_{1}^{\mathrm{ZAMS}}$, taken at
regular intervals in mass ratio $q$ between 0.1 and 1 assuming a flat
distribution (slope $\kappa=0$, e.g., Kouwenhoven et al., 2005). We consider
in Sec. 6 also variations with $\kappa=\pm 1$.

Finally, for each pair of masses
$(M_{1}^{\mathrm{ZAMS}},M_{2}^{\mathrm{ZAMS}})$ we choose $N_{p}$ different
periods $P^{\mathrm{ZAMS}}$ equally spaced in
$\log_{10}(P^{\mathrm{ZAMS}}/\mathrm{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_{1}<15\,M_{\odot}$ we assume a flat
distribution in $\log_{10}(P^{\mathrm{ZAMS}}/\mathrm{days})$ (Öpik, 1924;
Kobulnicky & Fryer, 2007); while for $M_{1}\geq 15\,M_{\odot}$ we assume a
powerlaw distribution in $\log_{10}(P^{\mathrm{ZAMS}}/\mathrm{days})$ with
exponent $\pi=-0.55$ (Sana et al., 2012). We explore values of $\pi=0\,$and
$\pi=-1$ (for all values of $M_{1}^{\mathrm{ZAMS}}$) 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 Sec. 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_{M_{1}}\times N_{q}\times N_{p}=50\times 50\times 100=250\,000$ binary
systems222This number is then multiplied by the number of natal kicks that we
draw for each core-collapse, see Sec. 2.2.1.. Increasing this resolution to
$N_{M_{1}}\times N_{q}\times N_{p}=100\times 100\times 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” 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
$\beta_{\mathrm{RLOF}}=\beta_{\mathrm{th}}$. Larger rates are likely to drive
the accretor out of thermal equilibrium333We do not model the internal
structure of the stars, therefore, we do not follow its possible bloating
during mass transfer, which could potentially enhance its mass loss and
related spin down. and lead to unstable mass transfer (Neo et al., 1977;
Hurley et al., 2002). For a large portion of the parameter space,
$\beta_{\mathrm{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 $\beta=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=\gamma_{\mathrm{RLOF}}J_{\mathrm{orb}}/(M_{\mathrm{don}}+M_{\mathrm{acc}})$,
with $\gamma_{\mathrm{RLOF}}=M_{\mathrm{don}}/M_{\mathrm{acc}}$, where
$J_{\mathrm{orb}}$ is the orbital angular momentum, $M_{\mathrm{don}}$ and
$M_{\mathrm{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.,
$\gamma_{\mathrm{RLOF}}=\sqrt{2}(M_{\mathrm{don}}+M_{\mathrm{acc}})^{2}/(M_{\mathrm{acc}}M_{\mathrm{don}})\equiv\gamma_{\mathrm{disk}}$,
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. $\gamma=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_{\mathrm{acc}}/M_{\mathrm{don}}<q_{\mathrm{crit}}$), we assume that the
system enters a common envelope phase. The threshold value $q_{\mathrm{crit}}$
is uncertain and depends on the evolutionary stage of the donor. In our
fiducial simulation $q_{\mathrm{crit,A}}=0.65$ for a main sequence (MS) donor
(de Mink et al., 2007), $q_{\mathrm{crit,B}}=0.4$ for a Hertzsprung gap donor
(Hurley et al., 2002), and $q_{\mathrm{crit,RSG}}=0.25$ for core-helium
burning and red supergiant donors (Claeys et al., 2014). We also consider
values of $q_{\mathrm{crit,A}}=0.25,\,0.8$, $q_{\mathrm{crit,B}}=1,\,0.5,\,0$,
and $q_{\mathrm{crit,RSG}}=1.0$ in our model variations.

We treat common envelope evolution using the $\alpha_{\rm
CE}\lambda$-formalism (Webbink, 1984; Livio & Soker, 1988; de Kool, 1990;
Hurley et al., 2002). In all variations, we use an analytic fit to the
$\lambda_{g}$ values of Dewi & Tauris (2000) for the binding energy parameter
$\lambda$. 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 $\alpha_{\rm 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
$\alpha_{\rm CE}=0.1$, i.e., inefficient use of the orbital energy to eject
the common envelope, and $\alpha_{\rm CE}=10$, intended as indicative of a
fairly extreme case of additional energy input. We also test the combination
of $\alpha_{\mathrm{CE}}=0.1$ and $10$ with $q_{\mathrm{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.2.1 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.

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

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

  3. 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 $\mathbf{v}_{k}$, drawn from a distribution
of amplitude $v_{k}\equiv|\mathbf{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 (TT98 hereafter). To re-calculate the post-CC
orbital parameters, their algorithm assumes instantaneous loss of the ejecta,
because the ejecta velocity $v_{\mathrm{ej}}$ is much larger than the orbital
velocity $v_{\mathrm{orb}}$ ($v_{\mathrm{orb}}\lesssim 10^{2}\,\mathrm{km\
s^{-1}}\ll v_{\mathrm{ej}}\sim 10^{4}\,\mathrm{km\ s^{-1}}$). 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_{2}=3.5\,M_{\odot}$
star, which is the most massive companion they considered. We also assume an
efficiency of momentum transfer444$\eta=0.5$ in the algorithm from TT98. 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, $\lesssim 10\,\mathrm{km\ s^{-1}}$ (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
$\alpha=10\,\mathrm{degrees}$ 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-\alpha\leq
45\,\mathrm{degrees}$. The kick amplitude $|\mathbf{v}_{k}|\equiv v_{k}$ is
drawn from a Maxwellian distribution with one-dimensional root mean squared
dispersion $\sigma_{\mathrm{kick}}=265\,\mathrm{km\ s^{-1}}$ (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 $\sigma_{\mathrm{kick}}=0,\,300,\,1000\,\mathrm{km\
s^{-1}}$ 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}\rightarrow v_{k}(1-f_{b})$, were $f_{b}$ is the fallback fraction taken
from555We correct for a missing parenthesis in their equations, which can be
found by dimensional analysis (Belczynski et al., 2017). Eq. 16 in Fryer et
al. (2012), corresponding to their “rapid SN engine” which reproduces the NS-
BH mass gap between $\sim 2-5\,M_{\odot}$ (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_{\odot}$.

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; Repetto & Nelemans 2015; Repetto et al. 2017, 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
$\sigma_{\mathrm{kick}}=100\,\mathrm{km\ s^{-1}}$ and no fallback down scaling
of the kick amplitude for the CC resulting in BH formation, and our fiducial
kick ($\sigma_{\mathrm{kick}}=265\,\mathrm{km\ s^{-1}}$, 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_{\odot}$
(Schwab et al., 2010; Knigge et al., 2011) drawing the natal kicks for these
NSs from a Maxwellian with $\sigma_{\mathrm{kick}}=30\,\mathrm{km\ s^{-1}}$
creating a double peaked distribution. For more massive NSs and BHs the kick
is drawn from a Maxwellian with $\sigma_{\mathrm{kick}}=265\,\mathrm{km\
s^{-1}}$.

In all our simulations, we draw $N_{\mathrm{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_{\mathrm{kick}}$ to 50.

###  2.3 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 Sec. 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 ($\sim$3 Myr, Zapartas et al., 2017a), 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

Figure 2: Evolution of an example binary system with
$M_{1}^{\mathrm{ZAMS}}=20\,M_{\odot}$, $M_{2}^{\mathrm{ZAMS}}=15\,M_{\odot}$
with $P^{\mathrm{ZAMS}}=100\,\mathrm{days}$, described in Sec. 3 (see also the
cartoon in Fig. 1). The top panel shows the primary, secondary, and total
mass. The bottom panel shows the 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.

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 into our
simulations. As a representative example, we consider a
$M_{1}^{\mathrm{ZAMS}}=20\,M_{\odot}$ primary star, a
$M_{2}^{\mathrm{ZAMS}}=15\,M_{\odot}$ secondary with $P^{\mathrm{ZAMS}}=100$
days (corresponding to a separation $a_{\mathrm{ZAMS}}\simeq 300\,R_{\odot}$).
Fig. 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.

###  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 $\sim$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_{\mathrm{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 $\Delta M_{2}\simeq 14\,M_{\odot}$, and becomes the more massive
star, now with a mass of about $28\,M_{\odot}$ (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 $\sim$$5.3\,M_{\odot}$ 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_{\odot}$ and a pre-SN orbital velocity of about
$v_{2}=(M_{1}/(M_{1}+M_{2}))v_{\mathrm{orb}}\simeq 12.5\,\mathrm{km\ s^{-1}}$.

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^{\mathrm{ZAMS}}=7$ days would experience a
very similar evolutionary path, but it would only widen to a pre-CC period of
$\sim$$45$ days, resulting in a pre-CC orbital velocity of the secondary
slightly larger than $30\,\mathrm{km\ s^{-1}}$.

Because of the rejuvenating effect of mass accretion, the secondary remains on
the MS for another $\sim$$3.6\,\mathrm{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 type666We do not attempt to distinguish between stripped SN types (IIb, Ib
or Ic depending on the presence or lack of He lines in the spectrum).
IIb/Ib/Ic (Filippenko, 1997), which forms a NS of about
$M_{\mathrm{NS}}=1.6\,M_{\odot}$, with a fallback fraction $f_{b}=0.13$. The
SN ejecta mass is $\Delta M_{\mathrm{SN}}=M_{1}^{\mathrm{pre-
CC}}-M_{\mathrm{NS}}\simeq 3.7\,M_{\odot}$, where $M_{1}^{\mathrm{pre-CC}}$ is
the primary mass just before the CC, and $M_{\mathrm{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, $M_{1}^{\mathrm{pre-CC}}+M_{2}^{\mathrm{pre-CC}}\simeq
33\,M_{\odot}$, 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_{\mathrm{k}}\gg
v_{\mathrm{orb}}$ results in disruption of the binary regardless of the
direction, (ii) an intermediate kick $v_{\mathrm{k}}\simeq v_{\mathrm{orb}}$
disrupts the system only if the kick orientation is favorable, while (iii) a
small kick, $v_{\mathrm{k}}\ll v_{\mathrm{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}\simeq 870,230$, and $43\,\mathrm{km\ s^{-1}}$ (respectively
corresponding to $1000,265$, and $50\,\mathrm{km\ s^{-1}}$ before the
rescaling of the kick amplitude with the fallback fraction). For comparison,
the pre-CC orbital velocity in our example system is
$v_{\mathrm{orb}}^{\mathrm{pre-CC}}\simeq 100\,\mathrm{km\ s^{-1}}$
(corresponding to $v_{2}^{\mathrm{pre-CC}}\simeq 12.5\,\mathrm{km\ s^{-1}}$).
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 Sec. 7 for a brief discussion). In the case of
disruption, the ejected star acquires a velocity $v_{\mathrm{dis}}\simeq
v_{2}^{\mathrm{pre-SN}}=12.5\,\mathrm{km\ s^{-1}}$. With
$v_{\mathrm{dis}}<30\,\mathrm{km\ s^{-1}}$, 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\,\mathrm{km\
s^{-1}}$), the angle $\theta$ between the kick direction and the pre-CC
orbital velocity of the primary, given by
$\cos\theta=\mathbf{v_{k}}\cdot\mathbf{v}_{1}^{\mathrm{pre-SN}}/(v_{k}\
v_{1}^{\mathrm{pre-SN}})$, determines whether the kick disrupts the system.
Kicks that are roughly aligned with the pre-CC orbital velocity, i.e.
$0\lesssim\theta\lesssim\pi/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
$v_{\mathrm{dis}}\simeq v_{2}^{\mathrm{pre-SN}}$.

With small kick amplitudes (here $v_{k}\lesssim 50\,\mathrm{km\ s^{-1}}$), 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 Sec. 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

Figure 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 $v_{2}$ (right panel), assuming
$q^{\mathrm{ZAMS}}=0.75$ and $\mu_{\mathrm{env}}=0.75$, which are
representative for the example system of Sec. 3. Blue colors correspond to a
decrease, red to an increase. The mass transfer efficiency
$\beta_{\mathrm{RLOF}}$ and the angular momentum parameter
$\gamma_{\mathrm{RLOF}}$ are assumed constant throughout the evolution.

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
$v_{\mathrm{dis}}\simeq v_{2}^{\mathrm{pre-CC}}$ (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_{2}$ of the
secondary star in the frame of the center of mass at any time during the
binary evolution as

$v_{2}=\frac{M_{1}}{M_{1}+M_{2}}v_{\mathrm{orb}}\equiv\frac{M_{1}}{M_{1}+M_{2}}\sqrt{\frac{G(M_{1}+M_{2})}{a}}\
\ ,$ (1)

where $G$ is the gravitational constant, $a$ is the semimajor axis of the
orbit (assumed to be circular), and $M_{1}$, $M_{2}$ are the two masses.

If we assume that a constant fraction
$\beta_{\mathrm{RLOF}}=|\dot{M}_{\mathrm{acc}}|/|\dot{M}_{\mathrm{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
$\gamma_{\mathrm{RLOF}}$ of the total specific orbital angular momentum, we
can relate the final and initial separations with

$\frac{a_{\mathrm{pre-
CC}}}{a_{\mathrm{ZAMS}}}=\left(\frac{M_{1}^{\mathrm{ZAMS}}}{M_{1}^{\mathrm{pre-
CC}}}\frac{M_{2}^{\mathrm{ZAMS}}}{M_{2}^{\mathrm{pre-
CC}}}\right)^{2}\left(\frac{M_{1}^{\mathrm{pre-CC}}+M_{2}^{\mathrm{pre-
CC}}}{M_{1}^{\mathrm{ZAMS}}+M_{2}^{\mathrm{ZAMS}}}\right)^{2\gamma_{\mathrm{RLOF}}+1}\
\ .$ (2)

The parameter $\beta_{\mathrm{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; Yoon et al.,
2017; Götberg et al., 2018). We parametrize the fraction of the total mass in
the envelope as $\mu_{\mathrm{env}}=M_{\mathrm{env}}/M$, or in other words
$1-\mu_{\mathrm{env}}$ is the fraction of mass in the helium core. The
parameter $\mu_{\mathrm{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\lesssim\mu_{\mathrm{env}}\lesssim
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
$\mu_{\mathrm{env}}$. Using this parametrization, we can relate the final
masses of the two star in the binary to their ZAMS masses as

$M_{1}^{\mathrm{pre-CC}}\simeq
M_{1}^{\mathrm{ZAMS}}-\mu_{\mathrm{env}}M_{1}^{\mathrm{ZAMS}}\ \ ,\\\
M_{2}^{\mathrm{pre-CC}}\simeq
M_{2}^{\mathrm{ZAMS}}+\beta_{\mathrm{RLOF}}\mu_{\mathrm{env}}M_{1}^{\mathrm{ZAMS}}\
\ .\hfill$ (3)

We can now express the final orbital velocity of the secondary
$v_{2}^{\mathrm{pre-CC}}$ in terms of the initial $v_{2}^{\mathrm{ZAMS}}$.
Introducing $q\equiv
q^{\mathrm{ZAMS}}=M_{2}^{\mathrm{ZAMS}}/M_{1}^{\mathrm{ZAMS}}$,
$\beta\equiv\beta_{\mathrm{RLOF}}$, $\mu\equiv\mu_{\mathrm{env}}$, and
$\gamma\equiv\gamma_{\mathrm{RLOF}}$ to simplify the notation, we can write

$\frac{v_{2}^{\mathrm{pre-
CC}}}{v_{2}^{\mathrm{ZAMS}}}=\left(\frac{q+\beta\mu-\mu
q-\beta\mu^{2}}{q}\right)^{2}\left(\frac{1+q}{q+1+\beta\mu-\mu}\right)^{\gamma+1}\
\ .$ (4)

In Fig. 3 we visualize the outcome of these analytic estimates for different
assumptions for the mass transfer efficiency
$\beta\equiv\beta_{\mathrm{RLOF}}$ and the angular momentum loss parameter
$\gamma\equiv\gamma_{\mathrm{RLOF}}$. To make this figure, we adopt the
parameters of the example system of Sec. 3 ($q=0.75,\
\mu_{\mathrm{env}}=0.75$).

The left panel of Fig. 3 shows the relative variation in the specific orbital
angular momentum
$h\stackrel{{\scriptstyle\mathrm{def}}}{{=}}J_{\mathrm{orb}}/(M_{1}+M_{2})=M_{1}M_{2}\sqrt{Ga}/(M_{1}+M_{2})^{3/2}$.
For fully conservative evolution ($\beta=1$, white region on the right of this
panel), $h$ stays constant. Also assuming $\gamma=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 $\gamma\lesssim 1$ result in a net increase of the specific orbital
angular momentum because mass is lost relatively faster than angular momentum,
while very large $\gamma\gtrsim 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 ($\gamma\gtrsim 2$) in non-conservative
($\beta\lesssim 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 $\beta$ and $\gamma$ 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
($\beta\lesssim 0.3$) experiencing large angular momentum losses
($\gamma\gtrsim 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, $\beta$ and
$\gamma$, 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 $\beta=1$. The expression for the final orbital velocity of the
secondary (Eq. 4) simplifies to

$\frac{v_{2}^{\mathrm{pre-
CC}}}{v_{2}^{\mathrm{ZAMS}}}=\left(\frac{(1-\mu)(q+\mu)}{q}\right)^{2}=\left(\frac{3+4q}{16q}\right)^{2}\
\ ,$ (5)

where we use as a typical value $\mu=3/4=0.75$ in the last step. Eq. 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.

$M_{1}^{\mathrm{ZAMS}}\geq 7.5\,M_{\odot}$
$q^{\mathrm{ZAMS}}\geq 0.1$
$0.15\leq\log_{10}(P^{\mathrm{ZAMS}}/\mathrm{days})\leq 5.5$$22_{-9}^{+26}\%$
$78_{-22}^{+9}\%$$1-\mathcal{D}=14_{-10}^{+22}\%$
7$36_{-15}^{+42}\%$$64_{-42}^{+15}\%$$\mathcal{D}=86_{-22}^{+10}\%$$25_{-10}^{+8}\%$
5.3$53_{-46}^{+12}\%$$23_{-20}^{+38}\%$$25_{-6}^{+23}\%$$75_{-8}^{+10}\%$
5.1$5_{-2}^{+14}\%$
$v_{\mathrm{dis}}\geq 30\,\mathrm{km\ s^{-1}}$$95_{-14}^{+2}\%$
$v_{\mathrm{dis}}<30\,\mathrm{km\ s^{-1}}$ Figure 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
$\mathcal{D}$. The errors on each fraction exclude the run without SN kicks
($\sigma_{\mathrm{kick}}=0\,\mathrm{km\ s^{-1}}$), which produces an
unrealistically low disruption fraction (cf. Table 1 and Sec. 6).

Using the typical values of our example system discussed in Sec. 3 ($q=0.75$),
we find $v_{2}^{\mathrm{pre-CC}}\simeq 0.25\,v_{2}^{\mathrm{ZAMS}}\simeq
20\,\mathrm{km\ s^{-1}}$. 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 $v_{2}^{\mathrm{pre-CC}}\simeq
22\,\mathrm{km\ s^{-1}}$, 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 Sec. 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.
Sec. 2). We describe the robustness of our findings against model assumptions
in Sec. 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 $\sim$$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).

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
($\mathcal{D}=86\%$ in our fiducial run), or the newly formed compact object
remains bound to the companion star. In Sec. 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 Sec. 5.3.

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

###  5.1 Ejected main sequence companions

Figure 5: Velocity distribution of MS stars ejected from a binary system at
the time of the first CC. The top panel shows the 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_{\odot}$ (i.e., roughly those that might
experience CC), and only MS stars more massive than $15\,M_{\odot}$ (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 $v_{\mathrm{dis}}\geq
30\,\mathrm{km\ s^{-1}}$. See also Fig. 14 for a wider velocity range, Fig. 13
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.

The fraction of CC events occurring with a MS companion in our fiducial run is
$\sim$76% (and it ranges from $69-90\%$ in our parameter variations). The vast
majority of these binaries which host a MS companion is 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,
$M_{2}^{\mathrm{post-CC}}\geq 7.5\,M_{\odot}$, where $M_{2}^{\mathrm{post-
CC}}$ 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 $M_{2}^{\mathrm{post-CC}}\geq
15\,M_{\odot}$, corresponding roughly to O-type stars. In the rest of this
section, we quote the ratios for the population of massive stars (i.e.
$M_{2}^{\mathrm{post-CC}}\geq 7.5\,M_{\odot}$) unless otherwise stated, see
also Table 1 for more information.

Although we do produce runaways, the majority of the systems ejects a star
slower than $30\,\mathrm{km\ s^{-1}}$, 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 $\mathcal{R}$ of walkaway to runaway stars produced per CC event in a
binary. From our fiducial simulation, we obtain
$\mathcal{R}_{\mathrm{7.5}}\simeq 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 Sec. 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\,\mathrm{km\ s^{-1}}$: 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. 14 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\,\mathrm{km\ s^{-1}}$. Massive stars
significantly faster than this are likely produced by different mechanisms
(see Boubert et al. 2017a, and 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\,\mathrm{km\ s^{-1}}$, 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_{\mathrm{dis}}\simeq 6\,\mathrm{km\
s^{-1}}$. Ejected MS companions with such velocity typically originate from
systems with $q_{\mathrm{ZAMS}}\gtrsim 0.5$ and long initial periods
($\log_{10}(P^{\mathrm{ZAMS}}/\mathrm{days})\simeq 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.

Figure 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_{\odot}$ and $15\,M_{\odot}$,
respectively. The vertical dashed lines mark the mean values of the
distributions.

Among low and intermediate mass walkaways with $v_{\mathrm{dis}}\simeq
6\,\mathrm{km\ s^{-1}}$, 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_{\odot}$. 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_{\odot}$, and $91\%$ of those more massive than $15\,M_{\odot}$ 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 $M_{2}^{\mathrm{ZAMS}}\simeq
0.5\,M_{1}^{\mathrm{ZAMS}}$, and we also assume that mass transfer is unstable
and leads to mergers if the mass ratio is too extreme
($M_{\mathrm{acc}}/M_{\mathrm{don}}<q_{\mathrm{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 $\sim$13 for masses larger $7.5\,M_{\odot}$.

Accounting for the durations of each evolutionary phase, we can also quantify
the fraction of MS stars more massive than $15\,M_{\odot}$ that are runaway or
walkaway stars. In our fiducial simulation, we find a runaway fraction for
masses larger than $15\,M_{\odot}$ of about $\sim$$0.5\%$ and a walkaway
fraction of $\sim$$10\%$. We only simulate binaries with primaries more
massive than $7.5\,M_{\odot}$, 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., 2017a). We show in Fig. 13 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 7: Equatorial velocity at the time of ejection for massive
($M_{\mathrm{dis}}\geq 7.5\,M_{\odot}$) 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 $\sim$$500\,\mathrm{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.

Figure 6 shows the distribution of distances traveled by massive, unbound
stars, calculated by multiplying their velocity $v_{\mathrm{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 star777The duration of the helium core burning depends on the
helium core mass, and thus the overshooting parameter. The value of $10\%$ of
the MS duration is typical for overshooting values larger than what assumed in
Pols et al. (1998).. 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_{\mathrm{dis}}\geq 7.5\,M_{\odot}$, we find the mean distance
$\langle L\rangle$ 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
$\langle L_{\mathrm{run}}\rangle=584$ pc considering only runaways (faster
than $30\,\mathrm{km\ s^{-1}}$) more massive than $7.5\,M_{\odot}$, but it
remains $\langle L_{\mathrm{walk}}\rangle=103$ pc considering massive
walkaways only (slower than $30\,\mathrm{km\ s^{-1}}$). For comparison, the
typical size of OB-associations is of the order of tens of pc, while the
Galactic thin and thick disks have a vertical scale height of $\sim$300 and
$\sim$1500 pc, respectively (e.g., Gilmore & Reid, 1983).

###  5.2 Spin of the main sequence ejected companions

Figure 7 shows the equatorial rotational velocity
$v_{\mathrm{dis}}^{\mathrm{eq}}$ for the ejected MS companions with
$M_{\mathrm{dis}}\geq 7.5\,M_{\odot}$. 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.

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_{\mathrm{dis}}<30\,\mathrm{km\ s^{-1}}$) 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
$v_{\mathrm{dis}}^{\mathrm{eq}}\simeq 200$ and $500\,\mathrm{km\ s^{-1}}$. 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_{\mathrm{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: 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. The top panel shows the corresponding cumulative
distributions.

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., 2017a). We distinguish three sub-populations: hydrogen-rich non-
degenerate secondaries (mostly of red-supergiant secondaries, corresponding to
$\sim$$4\%$ of the entire simulated population, including mergers and non-
disrupted binaries), white dwarfs (WD, $\sim$8% of the total), and helium
stars which have lost their hydrogen-rich envelope because of previous binary
interactions ($\sim 4\%$ of the total).

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 $\sim 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\lesssim v_{\mathrm{dis}}/\mathrm{[km\
s^{-1}]}\lesssim 90$, and helium stars are even faster with $70\lesssim
v_{\mathrm{dis}}/\mathrm{[km\ s^{-1}]}\lesssim 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_{\mathrm{dis}}<30\,\mathrm{km\ s^{-1}}$) helium stars. These are very
massive ($M_{2}\gtrsim 30\,M_{\odot}$ 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 Sec. 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
$\sigma_{\mathrm{kick}}=30\,\mathrm{km\ s^{-1}}$ for NSs less massive than
$1.35\,M_{\odot}$ (and our fiducial kick for more massive remnants), the
disruption fraction decreases to $\mathcal{D}=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
($\sigma_{\mathrm{kick}}=0\,\mathrm{km\ s^{-1}}$), then only
$\mathcal{D}=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 $f_{15}^{\mathrm{RW}}$ for masses larger than
$15\,M_{\odot}$ is low for all our parameter variations and never exceeds a
few percent. Similarly, the walkaway fraction $f_{15}^{\mathrm{WA}}$ 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 Sec. 6.3). Another case producing
slightly smaller $\mathcal{R}$ is fully non-conservative mass transfer
($\beta_{\mathrm{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
($\alpha_{\mathrm{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_{\mathrm{ZAMS}}\gtrsim 15\,M_{\odot}$, and the $q_{\mathrm{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 massive888The inspiraling star is assumed
not to accrete significantly during a common envelope event, which prevents
its growth by mass. runaway relatively small. Allowing for less stable late
RLOF (higher $q_{\mathrm{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 $\alpha_{\mathrm{CE}}$) are not sufficient to reconcile our
predictions with commonly-accepted inferences from observations (see Sec.
8.2).

The mean value of the mass of the walkaways, runaways, and all ejected
secondaries ($\langle M_{2}^{\mathrm{walk}}\rangle$, $\langle
M_{2}^{\mathrm{run}}\rangle$, and $\langle M_{2}\rangle$, respectively)
exceeds $7.5\,M_{\odot}$ in almost all our parameter variations. With fully
non-conservative mass transfer ($\beta_{\mathrm{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_{\mathrm{crit,B}}=1$), or if we take a negative slope of the mass ratio
distribution ($\kappa=-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 ($\sigma_{\mathrm{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_{\mathrm{crit,B}}$) impact
the velocity and mass distribution of the ejected stars. The sequence from
more stable to more unstable $q_{\mathrm{crit,B}}=1,\,0.5,\,0.4$ (fiducial
value), 0 shows progressively higher $\mathcal{R}$ indicating relatively fewer
runaways, but higher average velocity $\langle v\rangle$ (and consequently
also a larger $\langle L\rangle$). The increase in velocity and decrease in
average mass of the ejected companions for higher $q_{\mathrm{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 $\langle L\rangle\simeq 120$ pc, which corresponds to an average
distance for the runaways of $\langle L\rangle_{\mathrm{run}}\gtrsim 500$ pc,
while the walkaways alone reach $\langle L\rangle_{\mathrm{run}}\gtrsim 100$
pc.

These distances are significantly shortened only assuming fully conservative
mass transfer ($\beta_{\mathrm{RLOF}}=1$), or assuming higher angular momentum
loss (e.g., through a circumbinary disk,
$\gamma_{\mathrm{RLOF}}=\gamma_{\mathrm{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 $\langle L\rangle_{\mathrm{walk}}\simeq
66\,\mathrm{pc}$ is that with zero natal kicks
($\sigma_{\mathrm{kick}}=0\,{\mathrm{km\ s^{-1}}}$). 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_{\mathrm{dis}}\geq 7.5\,M_{\odot}$) 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

Figure 9: The mass function of massive ($M_{\mathrm{dis}}\geq
7.5\,M_{\odot}$), runaway (faster than $\geq 30\,\mathrm{km\ s^{-1}}$) stars
depends on the assumptions for the BH kick. The bottom panel shows our
fiducial simulation with $\sigma_{\mathrm{kick}}=265\,\mathrm{km\ s^{-1}}$,
including the fallback downscaling of the kick both for NSs and BHs. The
fallback is an important effect only for BHs. The central panel shows the mass
function for intermediate BH kicks ($\sigma_{\mathrm{kick}}=100\,\mathrm{km\
s^{-1}}$, no fallback downscaling) and our fiducial kick for NSs. The top
panel shows the distribution with our fiducial kick amplitude, but without
fallback downscaling (i.e., BH velocity kick). The dashed lines indicate the
$90^{\mathrm{th}}$ percentiles of the mass distributions.

Figure 10: Mass-velocity correlation for the three BH kick variations
discussed in Sec. 6.2. The brighter colors indicate where the probability per
system is higher on this plane, the colors have a logarithmic scale. The left
panel shows our fiducial distribution, the central panel shows the case for
intermediate BH kicks ($\sigma_{\mathrm{kick}}=100\,\mathrm{km\ s^{-1}}$ and
no fallback), and the right panel shows 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.

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_{\mathrm{dis}}\geq 7.5\,M_{\odot}$)
and fast runaways ($v_{\mathrm{dis}}\geq 30\,\mathrm{km\ s^{-1}}$) 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).

Fig. 9 shows the mass distribution of runaways
($v_{\mathrm{dis}}>30\,\mathrm{km\ s^{-1}}$) more massive than
$7.5\,M_{\odot}$ 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
$90^{\mathrm{th}}$ percentile of the massive runaway mass function is
$\sim$$25\,M_{\odot}$. 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.

The central panel of Fig. 9 shows the mass distribution assuming no fallback
downscaling, but a smaller kick amplitude distribution characterized by
$\sigma_{\mathrm{kick}}=100\,\mathrm{km\ s^{-1}}$ 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 $90^{\mathrm{th}}$ percentile of the mass distribution
shifts to $\sim$$30\,M_{\odot}$.

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
$90^{\mathrm{th}}$ percentile of the mass distribution is as high as
$\sim$$50\,M_{\odot}$.

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

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_{\mathrm{dis}}\gtrsim 70\,M_{\odot}$. 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 ($\sim 3\,\mathrm{Myr}$) above $\sim$$50\,M_{\odot}$. 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
$\sim$$70\,M_{\odot}$ 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\simeq Z_{\odot}$ (fiducial run), to
$Z=0.008$, $Z=0.0047$ to $Z=0.0002\simeq 10^{-2}Z_{\odot}$ a clear trend
emerges from Table 1. At lower metallicity, the fraction $\mathcal{D}$ 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 $\langle v\rangle$, which
goes from $11.2\,\mathrm{km\ s^{-1}}$ to $21.6\,\mathrm{km\ s^{-1}}$, but also
from the progressive decrease of the ratio of walkaways per runaway
$\mathcal{R}$, regardless of the lower mass cut. Nevertheless, walkaways
outnumber runaways at all the metallicities we explore.

In our lowest metallicity run, the runaways ($v_{\mathrm{dis}}\geq
30\,{\mathrm{km\ s^{-1}}}$) 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.

Table 1: Outcome of our model variations. Average distance traveled from the ejection to CC for all massive companions ejected, for massive runaways only, and massive walkaways only ($\langle L\rangle$, $\langle L_{\mathrm{run}}\rangle$, $\langle L_{\mathrm{walk}}\rangle$), average and median velocity $\langle v\rangle$ and $v_{\mathrm{med}}$, average masses for ejected main sequence companions, runaways only, and walkaways only ($\langle M_{2}\rangle$, $\langle M_{2}^{\mathrm{run}}\rangle$, $\langle M_{2}^{\mathrm{walk}}\rangle$), produced ratio $\mathcal{R}$ of walkaways to runaways for all masses, above $7.5\,M_{\odot}$, and above $15\,M_{\odot}$, fraction of all the binaries simulated resulting in a merger $\mathcal{M}$, disruption fraction among the non-merging binaries $\mathcal{D}$, and fraction of runaways and walkaways with masses larger than $15\,M_{\odot}$ from the disruption of binaries in a steady state population ($f^{\mathrm{RW}}_{15}$, $f^{WA}_{15}$).  Physical Assumptions | Parameter | value | $\langle L\rangle$ | $\langle L_{\mathrm{run}}\rangle$ | $\langle L_{\mathrm{walk}}\rangle$ | $\langle v\rangle$ | $v_{\mathrm{med}}$ | $\langle M_{2}\rangle$ | $\langle M_{2}^{\mathrm{run}}\rangle$ | $\langle M_{2}^{\mathrm{walk}}\rangle$ | $\mathcal{R}_{\mathrm{MS}}$ | $\mathcal{R}_{7.5}$ | $\mathcal{R}_{15}$ | $\mathcal{M}$ | $\mathcal{D}$ | $f^{\mathrm{RW}}_{\mathrm{15}}$ | $f^{\mathrm{WA}}_{\mathrm{15}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
$\mathrm{[pc]}$ | $\mathrm{[pc]}$ | $\mathrm{[pc]}$ | $\mathrm{[km\ s^{-1}]}$ | $\mathrm{[km\ s^{-1}]}$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | [%] | [%]
Fiducial population |  see Sec. 2 | 126 | 584 | 103 | 12.4 | 10.4 | 10.8 |  09.8 | 10.9 | 20.1 | 19.9 | 26.7 | 0.22 | 0.86 | 0.5 | 10.1
Mass transfer efficiency | $\beta_{\mathrm{RLOF}}$ | 0 | 139 | 493 |  098 | 15.1 | 11.6 | 7.7 | 8.6 |  07.6 | 10.3 |  08.7 |  05.8 | 0.23 | 0.86 | 0.3 |  01.5
0.5 | 122 | 363 | 111 | 12.8 | 10.5 | 10.2 | 12.0 | 10.1 | 21.7 | 20.7 | 11.0 | 0.22 | 0.87 | 1.2 |  08.6
1 |  082 | 214 |  080 | 11.2 |  09.8 | 12.1 | 11.9 | 12.1 | 33.0 | 46.9 | 26.6 | 0.21 | 0.87 | 0.7 | 14.7
Angular momentum loss | $\gamma_{\mathrm{RLOF}}$ | $\gamma_{\mathrm{disk}}$ |  072 | 544 |  064 | 10.7 |  09.8 | 10.5 |  06.4 | 10.7 | 22.6 | 59.0 | 53.1 | 0.36 | 0.85 | 0.2 |  07.3
1 | 129 | 592 | 104 | 12.6 | 10.5 | 10.8 | 10.0 | 10.9 | 19.3 | 18.5 | 25.3 | 0.22 | 0.86 | 0.6 |  09.9
Common envelope efficiency | $\alpha_{\mathrm{CE}}$ | 0.1 | 126 | 581 | 103 | 12.4 | 10.3 | 10.9 | 11.3 | 10.9 | 23.8 | 19.8 | 26.3 | 0.23 | 0.86 | 0.5 | 10.1
10 | 135 | 730 | 103 | 12.7 | 10.5 | 10.6 |  07.5 | 10.9 | 12.6 | 18.7 | 26.5 | 0.16 | 0.84 | 0.5 | 10.0
Mass ratio for case A merger | $q_{\mathrm{crit,\ A}}$ | 0.80 | 128 | 593 | 104 | 12.3 | 10.3 | 10.7 |  09.6 | 10.7 | 20.1 | 19.9 | 27.5 | 0.24 | 0.86 | 0.5 | 10.2
0.25 | 134 | 637 | 102 | 13.0 |  09.4 | 11.1 | 10.3 | 11.1 | 17.5 | 15.9 | 23.0 | 0.17 | 0.86 | 0.6 |  09.4
Mass ratio for case B merger | $q_{\mathrm{crit,\ B}}$ | 1.0 | 100 | 407 | 100 |  09.5 | 10.8 |  08.6 |  03.5 |  08.8 | – | – | – | 0.43 | 0.89 | – |  05.0
0.5 | 103 | 717 | 102 |  09.1 |  09.3 |  08.3 |  03.0 |  08.4 | – | – | – | 0.48 | 0.89 | – |  05.0
0.0 | 167 | 870 | 107 | 13.8 | 10.1 | 10.7 | 9.4 | 10.8 | 12.4 | 11.6 | 24.5 | 0.16 | 0.85 | 0.6 | 10.1
Mass ratio for case C merger | $q_{\mathrm{crit,\ RSG}}$ | 1.0 | 123 | 509 | 82 | 14.6 | 11.4 | 9.9 | 8.1 | 10.2 | 7.9 | 9.5 | 17.0 | 0.31 | 0.80 | 0.6 | 6.3
Mass ratio for case C merger and | $(q_{\mathrm{crit,\ RSG}},\alpha_{\mathrm{CE}})$ | (1.0, 0.1) | 109 | 626 | 79 | 13.3 | 10.5 | 10.4 | 12.3 | 10.3 | 23.7 | 17.7 | 19.4 | 0.41 | 0.82 | 0.5 | 6.4
Common Envelope efficiency | (1.0, 10) | 177 | 764 | 91 | 17.9 | 13.2 | 9.6 | 7.9 | 9.9 | 6.0 | 6.8 | 12.2 | 0.18 | 0.78 | 0.7 | 6.4
Natal kick amplitude | $\sigma_{\mathrm{kick}}\ [\mathrm{km\ s^{-1}}]$ | 0 |  066 | – |  066 |  03.8 | 10.2 |  01.6 | – |  01.6 | – | – | – | 0.22 | 0.16 | – | –
300 | 128 | 572 | 103 | 12.7 | 10.4 | 11.0 | 10.3 | 11.0 | 18.7 | 18.0 | 22.9 | 0.22 | 0.87 | 0.6 | 10.3
1000 | 132 | 490 | 102 | 13.9 | 10.9 | 11.7 | 13.8 | 11.5 | 13.8 | 12.0 | 11.6 | 0.22 | 0.91 | 1.2 | 11.2
BH kicks $\sigma_{\mathrm{kick}}=\left\\{\phantom{\frac{0}{0}}\right.$ |  $100\,\mathrm{km\ s^{-1}},\ (f_{b}=0$) | for BH | 121 | 395 | 97 | 13.9 | 10.7 | 12.8 | 18.7 | 12.4 | 13.7 | 11.3 | 9.0 | 0.22 | 0.97 | 0.1 | 12.1
$265\,\mathrm{km\ s^{-1}}$ | for NS
Double maxwellian with $\sigma_{\mathrm{kick}}=\left\\{\phantom{\frac{0}{0}}\right.$ | $30\,\mathrm{km\ s^{-1}}$ |  for $M_{\mathrm{NS}}\leq 1.35\,M_{\odot}$ | 105 | 423 | 087 | 12.2 | 09.5 | 10.2 | 14.5 | 10.0 | 30.0 | 17.5 | 15.5 | 0.22 | 0.65 | 0.5 | 04.9
$265\,\mathrm{km\ s^{-1}}$ |  for $M_{\mathrm{NS}}>1.35\,M_{\odot}$
Restricted kick directions | $\alpha<10\,\mathrm{deg}$ | 127 | 582 | 103 | 12.5 | 10.5 | 10.9 |  09.7 | 11.0 | 19.8 | 19.4 | 27.4 | 0.22 | 0.87 | 0.6 | 10.3
$90\,\mathrm{deg}-\alpha<45\,\mathrm{deg}$ | 128 | 572 | 102 | 12.5 | 10.4 | 10.8 | 10.7 | 10.8 | 18.1 | 17.2 | 19.4 | 0.22 | 0.86 | 0.5 | 10.0
Fallback fraction | $f_{b}$ | 0 | 121 | 395 |  097 | 13.9 | 10.7 | 12.8 | 18.7 | 12.4 | 13.7 | 11.3 | 9.0 | 0.22 | 0.97 | 1.5 | 12.1
Period distribution slope | $\pi$ | -1 | 128 | 581 | 102 | 12.9 | 10.7 | 10.6 | 10.1 | 10.6 | 18.7 | 17.6 | 22.8 | 0.28 | 0.89 | 0.6 |  08.9
0 | 120 | 593 | 102 | 11.5 |  09.8 | 10.8 |  09.0 | 10.9 | 23.8 | 26.4 | 35.9 | 0.17 | 0.85 | 0.4 | 10.7
Initial mass function slope | $\alpha$ | -1.9 | 125 | 554 | 101 | 12.8 | 10.5 | 11.8 | 10.9 | 11.8 | 18.9 | 18.0 | 23.9 | 0.24 | 0.80 | 0.6 |  09.1
-3 | 125 | 626 | 104 | 11.7 | 10.1 |  09.6 |  08.3 |  09.6 | 22.6 | 24.2 | 32.7 | 0.20 | 0.93 | 0.5 | 11.6
Mass ratio slope | $\kappa$ | -1 | 163 | 629 | 130 | 13.2 | 11.1 |  07.8 |  06.8 |  07.9 | 14.5 | 14.5 | 20.7 | 0.36 | 0.86 | 0.4 |  06.3
1 |  099 | 526 |  084 | 11.7 |  09.7 | 13.0 | 12.7 | 13.1 | 28.3 | 27.4 | 33.2 | 0.14 | 0.87 | 0.5 | 12.4
Metallicity | $Z$ | 0.0002 | 279 | 706 | 163 | 21.6 | 16.9 | 9.2 | 10.8 |  08.9 |  04.4 |  03.7 | 3.0 | 0.15 | 0.77 | 2.6 |  07.7
0.0047 | 190 | 718 | 130 | 15.9 | 12.5 | 10.5 | 11.0 | 10.4 | 10.7 | 8.8 | 11.1 | 0.19 | 0.84 | 1.2 | 10.3
0.008 | 164 | 682 | 123 | 14.4 | 11.9 | 10.5 |  09.9 | 10.6 | 13.7 | 12.8 | 17.4 | 0.21 | 0.84 | 0.1 | 11.0
0.03 | 109 | 576 |  087 | 11.6 |  09.7 | 11.2 |  09.9 | 11.3 | 22.5 | 21.6 | 32.6 | 0.24 | 0.88 | 0.5 | 10.0
Initial spin distribution | R-A15 | 126 | 582 | 103 | 12.4 | 10.4 | 10.9 |  09.8 | 10.9 | 20.1 | 19.8 | 26.5 | 0.22 | 0.86 | 0.5 | 10.1
Figure 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. The left panel shows the distribution for
NS companions, the right panel shows 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. The top panels show the corresponding cumulative distributions.

##  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
$1-\mathcal{D}=14_{-10}^{+14}\%$ of the total (see Sec. 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 $\sim$$66\,\mathrm{km\ s^{-1}}$, much lower than
the average of $\sim$$330\,\mathrm{km\ s^{-1}}$ (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_{\mathrm{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_{\mathrm{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.

For NS+MS binaries our fiducial simulation shows a median systemic velocity of
$\sim 20{\mathrm{km\ s^{-1}}}$, 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-\mathcal{D}\simeq 35$%. It
primarily adds systems with low systemic velocities. The median of the
normalized distribution for this simulation lies near $7{\mathrm{km\
s^{-1}}}$.

In our fiducial simulations for BH+MS binaries we obtain a bimodal
distribution with a large peak at $0{\mathrm{km\ s^{-1}}}$ accounting for
nearly 80% of all systems and a second broader component peaking at $\sim
10{\mathrm{km\ s^{-1}}}$. 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 $\sim$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_{\mathrm{sys}}\simeq 0\,{\mathrm{km\ s^{-1}}}$ 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-\mathcal{D}\simeq 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{\mathrm{km\
s^{-1}}}$.

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 $M_{2}^{\mathrm{pre-CC}}>M_{1}^{\mathrm{pre-CC}}$) because of the
previous mass transfer phase. In our fiducial run, for the systems remaining
bound the average pre-CC mass ratio is $\langle q_{\mathrm{pre-
CC}}\rangle\equiv\langle M_{2}^{\mathrm{pre-CC}}/M_{1}^{\mathrm{pre-
CC}}\rangle\simeq 2.9$. The majority of the bound systems, about $60\%$ for
all variations with non-zero natal kicks, have separations of less than about
$500\,R_{\odot}$. 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 $\sim 6000\,R_{\odot}$.

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^{-4}$). 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\,{\mathrm{km\ s^{-1}}}$ (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
$\sim$ $5\,{\mathrm{km\ s^{-1}}}$ 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_{\odot}$ (red) with a Gaussian distribution with full-width-half-
maximum of 10 ${\mathrm{km\ s^{-1}}}$, 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\,{\mathrm{km\ s^{-1}}}$, 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\,{\mathrm{km\ s^{-1}}}$, the effect is insufficient to reconcile the
observed runaway fraction with our results (see Sec. 8.2). The ratio of
walkaways to runaways more massive than $15\,M_{\odot}$ remains
$\mathcal{R}_{15}^{\mathrm{SFH}}\sim 15$ even accounting for the velocity
dispersion of the parent association (cf. $\mathcal{R}_{15}\sim 27$ before
accounting for it).

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 $\delta$-distribution with a velocity
$\mathbf{v}_{\mathrm{lag}}=(v_{\mathrm{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 $\delta$-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_{\mathrm{lag}}\simeq 5\,{\mathrm{km\ s^{-1}}}$ (Reid et al.,
2014) would barely change the fraction of _apparent_ O-type runaways we
predict. To infer an O-type runaway fraction of $\sim$10%, the true lag
velocity of high-mass star forming regions would need to be larger than
$\gtrsim$$15\,{\mathrm{km\ s^{-1}}}$, and would need not to have been
corrected when making inferences from observations.

We note that for O-type runaways with lifetimes of $\lesssim 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.

Figure 12: The yellow histogram shows the convolution of the distribution of
ejected MS stars more massive than $15\,M_{\odot}$ (red histogram, cf. Fig. 5)
with a Gaussian with full-width-half-maximum of $10\,{\mathrm{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.

###  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 $\sim$$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 preparation), 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:

$f^{\mathrm{RW}}\simeq
f_{\mathrm{bin}}\times(1-f_{\mathrm{mergers}})\times\mathcal{D}\times(1+f_{\mathrm{acc}})^{\alpha}\times
f_{v>30}\ \ ,$ (6)

where we can assume a binary fraction $f_{\mathrm{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_{\mathrm{mergers}}\simeq 0.8$ (Sana et al., 2012), the disruption
fraction is $\mathcal{D}\simeq 0.8$ (this study, but also, e.g., Eldridge et
al., 2011), the fraction of material accreted, for which we assume
$f_{\mathrm{acc}}\simeq 0.2$ (Packet, 1981), and the fraction of ejected
companions fast enough to be a runaway which we can estimate999We neglect here
that $\mathcal{D}$ is the disruption fraction including also binaries with no
MS stars at the pre-CC stage from $\mathcal{D}$ and $\mathcal{R}_{15}$, and is
$f_{v>30}\simeq\mathcal{D}/(1+\mathcal{R}_{15})\simeq 0.04$ for our fiducial
run. Assuming an IMF slope $\alpha=-2.3$ we obtain a runaway fraction of
$f^{\mathrm{RW}}\simeq 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. Sec. 8.1), or any combination
of the previous.

Since $v_{\mathrm{dis}}\simeq v_{2}^{\mathrm{pre-CC}}\propto
v_{\mathrm{orb}}^{\mathrm{pre-CC}}\propto J_{\mathrm{orb}}$, with
$J_{\mathrm{orb}}$ orbital angular momentum, to shift the peak of the velocity
distribution in Fig. 5 from $\sim$$6\,\mathrm{km\ s^{-1}}$ to
$v_{\mathrm{dis}}\gtrsim 30\,\mathrm{km\ s^{-1}}$, 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 $\lesssim$10 $\mathrm{km\
s^{-1}}$ (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_{\mathrm{dis}}\geq 7.5\,M_{\odot}$. 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 ($\lesssim 50\,\mathrm{Myr}$) stars within 3 Kpc of the
Solar neighborhood using two Maxwellians: one for the “low-velocity”
($v\lesssim 30\,\mathrm{km\ s^{-1}}$) and one for the “high-velocity” group.
Our findings suggest that $\sim$$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 $\tau_{\mathrm{SNR}}\simeq 10^{4}$ years, the typical
distance traveled by an ejected companion is of order $D\equiv
v_{\mathrm{dis}}\times\tau_{\mathrm{SNR}}\simeq 0.1\,\mathrm{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. (2017) 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 $\sim$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\,\mathrm{km\ s^{-1}}$ 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
$\sim$$400\,\mathrm{km\ s^{-1}}$ and $\sim$$150\,\mathrm{km\ s^{-1}}$ 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 $\mathrm{km\
s^{-1}}$ are only marginally reached in our simulations (cf. Fig. 5, Fig. 8,
and Fig. 14). 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 $\sim 70\,\mathrm{km\ s^{-1}}$, 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; Repetto et al., 2012; Fragos et al., 2013; 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 $\sim 300\,R_{\odot}$, and most bound system have
separations smaller than $\sim 500\,R_{\odot}$, 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. Sec. 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 $\mathrm{km\ s^{-1}}$, and since the secondary star is also an accretor and
rejuvenated, they are able to travel a distance comparable to $\langle
L_{\mathrm{walk}}\rangle$. 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_{\mathrm{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
Sec. 3 and Sec. 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 is disrupted at the first CC (with a few
notable exceptions, see also Sec. 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. 2016b), but it might possibly
decrease the efficiency of their mechanical feedback because of the lower