Core-Collapse Supernovae in Binaries as the Origin of Galactic Hyper-Runaway Stars

Several stars detected moving at velocities near to or exceeding the Galactic escape speed likely originated in the Milky Way disc. We quantitatively explore the `binary supernova scenario' hypothesis, wherein these stars are ejected at large peculiar velocities when their close, massive binary companions undergo a core-collapse supernova and the binary is disrupted. We perform an extensive suite of binary population synthesis simulations evolving massive systems to determine the assumptions and parameters which most impact the ejection rate of fast stars. In a simulation tailored to eject fast stars, we find hyper-runaway star progenitor binaries composed of a massive ($\sim$30 M$_{\odot}$) primary and a $\sim$3-4 M$_{\odot}$ companion on an orbital period that shrinks to $\lesssim$1 day prior to the core collapse following a common envelope phase. The black hole remnant formed from the primary must receive a natal kick $\gtrsim$1000 km s$^{-1}$ to disrupt the binary and eject the companion at a large velocity. We compare the fast stars produced in these simulations to a contemporary census of early-type Milky Way hyper-runaway star candidates. We find that these rare objects may be produced in sufficient number only when poorly-constrained binary evolution parameters related to the strength of post-core collapse black hole natal kicks and common envelope efficiency are adjusted to values currently unsupported -- but not excluded -- by the literature. We discuss observational implications that may constrain the existence of these putative progenitor systems.


INTRODUCTION
In recent years, works have reported with an increasing frequency detections of early-type main sequence (MS) stars in the Galactic halo moving at very high velocities, near to or exceeding the Galactic escape velocity at their position (e.g., Brown et al. 2005;Hirsch et al. 2005;Edelmann et al. 2005;Brown et al. 2006Brown et al. , 2009Brown et al. , 2012Brown et al. , 2014Zhong et al. 2014;Huang et al. 2017;Irrgang et al. 2019;Koposov et al. 2019, and references therein). For reference, the Galactic escape velocity in the Solar Neighbourhood is ∼530 km s −1 (Piffl et al. 2014;Williams et al. 2017) and falls to 400 km s −1 beyond 50 kpc from the centre of the Milky Way (Williams et al. 2017). The population of high velocity stars is increasing further in the Gaia era (Gaia Collaboration et al. 2016Collaboration et al. , 2018, with the European Space Agency satellite providing precise astrometry for billions of Milky Way sources. The E-mail: evans@strw.leidenuniv.nl short-lived nature of early-type stars and the dearth of star formation in the Galactic stellar halo suggest that these fast stars were not likely formed in-situ in the halo; rather, they were likely accelerated and ejected from their primal birthplaces. They therefore flag extreme astrophysical and dynamical processes occurring in the Milky Way. These processes as well as the uncertain initial conditions and physics governing them can be elucidated by studying the properties and kinematics of these fast-moving stars. See  for a recent review on these objects. There exist a number of mechanisms to accelerate stars to such extreme velocities. Hills (1988) predicted that tight stellar binaries in the Galactic Centre (GC) could be tidally disrupted by a supermassive black hole (SMBH) lurking in the centre of the Milky Way. One member of the binary is captured in orbit around the SMBH while its companion is ejected with a very high velocity, ∼1000 km s −1 , giving rise to a population of so-called 'hyper-velocity stars' (HVS). Variations on this mechanism include the interaction of a single star with a binary system consisting of two SMBHs or an SMBH and an intermediate mass black hole (e.g. Yu & Tremaine 2003;Gualandris et al. 2005;Sesana et al. 2006;Guillochon & Loeb 2015), or the disruption of a globular cluster by a single SMBH or binary massive black hole pair in the GC (Capuzzo-Dolcetta & Fragione 2015;Fragione & Capuzzo-Dolcetta 2016). Regardless, a GC origin for hypervelocity stars is a shared element of the above mechanisms. Theoretical estimates place the HVS ejection rate from the GC on the order of 10 −4 yr −1 (Hills 1988;Perets et al. 2007;Zhang et al. 2013).
With perhaps the notable exception of S5-HVS1 (Koposov et al. 2019), the GC cannot be identified indisputably as the origin of many HVS candidates with contemporary astrometric measurements. However, especially given the high-precision astrometry provided by the second data release (DR2) of the Gaia mission (Gaia Collaboration et al. 2016Collaboration et al. , 2018, the GC can be ruled out as the spatial origin of many HVS candidates (see e.g. Irrgang et al. 2018). Other mechanisms must therefore be invoked to explain the extreme velocity of these stars. While tidal disruption of infalling dwarf galaxies (Abadi et al. 2009) or ejection from the Large Magellanic Cloud (Przybilla et al. 2008a;Boubert & Evans 2016;Boubert et al. 2017;Lennon et al. 2017;Erkal et al. 2019) or M31 (Sherwin et al. 2008) can explain extreme velocity stars whose trajectories seem to point towards extragalactic space, the most plausible origin for many high velocity star candidates seems to be the Milky Way disc (e.g., Heber et al. 2008b;Silva & Napiwotzki 2011;Palladino et al. 2014;Irrgang et al. 2018Irrgang et al. , 2019Marchetti et al. 2019). This is in spite of the fact that theoretical studies predict GC-ejected high-velocity stars to far outnumber disc-ejected high-velocity stars (Bromley et al. 2009;Kenyon et al. 2014).
A number of processes can be invoked to explain the existence of these disc-ejected high velocity stars. In a tight white dwarf/helium star binary, the deposition of a critical amount of helium onto the accreting white dwarf can result in the thermonuclear detonation of the white dwarfa proposed progenitor for Type Ia supernovae (e.g. Justham et al. 2009). The donor star can be ejected at a velocity unbound to the Milky Way and therefore be observed as a hyper-velocity helium star or (eventually) a hyper-velocity white dwarf (Hansen 2003;Geier et al. 2013Geier et al. , 2015. For main sequence high-velocity stars seemingly ejected from the Galactic disc, on the other hand, two main processes are typically blamed. In the dynamical ejection scenario (DES, e.g. Poveda et al. 1967;Leonard & Duncan 1990;Leonard 1991;Perets &Šubr 2012;Oh & Kroupa 2016), exchange encounters in dense stellar systems (Aarseth 1974) may eject stars at high velocities. In the binary supernova scenario (BSS; e.g. Blaauw 1961;Boersma 1961;Tauris & Takens 1998;Portegies Zwart 2000;Tauris 2015;Renzo et al. 2019b), the massive primary in a tight binary explodes in a core-collapse (CC) supernova, disrupting the binary and ejecting its main-sequence companion with a velocity comparable to its pre-CC orbital velocity. Both processes are known to occur in the Milky Way (Hoogerwerf et al. 2001;Jilinski et al. 2010) and are generally thought to be responsible for the known sample of 'runaway stars' with v ≥ 30 − 40 km s −1 (Blaauw 1961), though their relative contribution is not yet well-constrained (see Hooger-werf et al. 2001;Renzo et al. 2019b). With characteristic ejection speeds on the order of a few tens of km s −1 , it is not yet known whether these mechanisms can eject stars on the order of hundreds of km s −1 with sufficient frequency to explain the current known sample of 'hyper-runaway stars' (HRSs) -runaway stars ejected near to or ab ove the Galactic escape velocity at their location. While ejection velocities in the neighbourhood of ∼1000 km s −1 are possible in both the DES (Leonard 1991) and BSS (e.g., Tauris & Takens 1998;Tauris 2015) scenarios, these situations are thought to be rare. Recent N-body simulations of young star clusters have found that ejections in excess 200 km s −1 are very rare (Perets &Šubr 2012;Oh & Kroupa 2016). Binary population synthesis models simulating a large number of binary systems show that ejections above 200 km s −1 are vastly outnumbered by ejections on the order of ∼10 km s −1 (Portegies Zwart 2000; Eldridge et al. 2011;Renzo et al. 2019b).
Extreme velocity stars are interesting beyond their status as astrophysical oddities -the violent and uncertain physical processes that generate them leave an imprint on their kinematics and properties. The population of stars ejected via the BSS and their mass and velocity distributions provide constraints on many uncertain parameters governing binary evolution and core-collapse supernovae -in particular the debated physics of the common envelope phase and the nature of the natal 'kicks' imparted on compact objects produced following CC events. Stars ejected via the DES can reveal information about the initial conditions describing their parent clusters. Genuine hyper-velocity stars ejected from the GC offer a convenient 'back door' into studying the dustobscured and source-crowded GC environment (e.g., Zhang et al. 2013;Madigan et al. 2014;Rossi et al. 2017), in particular the origin and nature of the Milky Way's nuclear star cluster (see Böker 2010, for a review), the interplay between Sgr A* and its environment (see Genzel et al. 2010, for a review), and the growth of Sgr A* via tidal disruption of former HVS companions (Bromley et al. 2012). The long journey from GC to outer halo makes GC-ejected hypervelocity stars intriguing dynamical tracers for studying the size, mass and shape of the Galactic dark matter halo (e.g., Gnedin et al. 2005;Yu & Madau 2007;Kenyon et al. 2008Kenyon et al. , 2014Rossi et al. 2017;Contigiani et al. 2019).
In this paper, we focus on HRSs ejected following a core-collapse event occurring in a binary system, building on the recent work by Renzo et al. (2019b), hereafter R19. They use a rapid binary population synthesis code to examine the kinematic signatures of ejected stars. To account for uncertainties in the initial conditions and physical processes important for binary evolution, they ran an extensive grid of simulations varying relevant assumptions. This allowed them to i) make concrete observational predictions by showing which aspects of the kinematic signatures are robust to model variations and ii) identify which assumptions and processes most strongly impact the kinematic observables, so that they may be constrained by future observations. They find in general across all model variations that the majority of binaries disrupted following the CC of the primary eject slow-moving 'walkaway stars' (v 30 km s −1 ), a term first introduced in de Mink et al. (2012). They also find, however, that the absolute number of unbound companions and their velocity distribution depend intimately on assumptions about the natal kick delivered to remnant neutron stars and black holes (BHs) from the asymmetric supernova explosion. While high-velocity ejections were not the main focus of the work, R19 remark on the presence of a minor peak in the ejection velocity distribution in the vicinity of 100 km s −1 v 400 km s −1 (c.f. Fig. C.2 in R19), resulting from binaries which experienced common envelope evolution without significant accretion by the secondary.
In this work we employ the same simulation framework and approach as R19 to focus further on these fastest ejections: we study whether the BSS is able to eject fast stars with sufficient frequency to match the current sample of identified HRS candidates in the Milky Way. We vary relevant uncertain physical parameters and determine which most strongly affect the likelihood of a system ejecting a companion at a very high velocity.
This paper is organized as follows. In Sec. 2 we describe our binary population synthesis code, the physical assumptions we use in all our models and how we compare the simulation results to current HRS candidates. In Sec. 3 we outline how we build and characterize a sample of known HRS candidates from the literature. In Sec. 4 we show how the number of observable HRSs predicted from our model runs differ from the observed sample. In Sec. 5 we discuss the implications and limitations of our findings, and our conclusions can be found in Sec. 6.

BINARY POPULATION SYNTHESIS MODELS
We generate and evolve mock populations of massive binary systems using the population synthesis and evolution code binary_c (Izzard et al. 2004(Izzard et al. , 2006(Izzard et al. , 2009, based on the binary star evolution (BSE) code of Hurley et al. (2002) and the algorithms of Tout et al. (1997). BSE itself uses the analytical fits of Hurley et al. (2000) to evolutionary models of single stars from Pols et al. (1998). The current version of binary_c also includes updates from de Mink et al. (2013), Schneider et al. (2014), and Claeys et al. (2014) for improved models of stellar rotation, stellar lifetimes, and Roche lobe overflow mass transfer rates respectively. We first generate a mock fiducial population, taking the same initial conditions and free physical parameters of R19's fiducial model 1 , based on reasonable assumptions consistent with the current theoretical and observational understanding of the field. We briefly explain these parameters in the next subsection and refer the reader to R19 for more information. Following the approach of R19, we vary relevant physical parameters one-by-one to examine their impact on the number of ejected high-velocity stars. As we will see, however, the number of high-velocity ejections in the fiducial simulation is quite small. While varying individual parameters can provide valuable hints towards which parameters are most important for high-velocity ejections, the impacts of varying individual parameters are difficult to verify robustly when comparing such small populations. We therefore instead run an 'optimized' simulation with all relevant free parameters optimistically tuned in directions which favour 1 This simulation is available at http://cdsarc.u-strasbg.fr/ ftp/J/A+A/624/A66/files/A_fiducial.sn.dat.gz.
high-velocity ejections of companions. By setting parameters back to their fiducial values one-by-one, we can more confidently assess their effect on the frequency of high-velocity ejections. See Table 1 for a summary of the initial conditions and physical parameters assumed in each simulation.

Fiducial Model Initial Conditions and Physical Assumptions
Here we briefly define and state our assumptions for a number of initial conditions and binary evolution parameters in the fiducial model. Each binary in our fiducial simulation is described by a zero-age main sequence (ZAMS) mass of the primary (M1), a ZAMS mass ratio between stars (q ≡ M2/M1) and an initial orbital period P . We assume that these parameters are independent from each other, although see Moe & Di Stefano (2017). We generate systems on a grid in this space and later statistically weigh each system, assigning each a probability according to our assumed distributions of these parameters, pi ∝ f (M1)f (q)f (P ) (see also R19, Appendix B). We select primary masses logarithmically-spaced in the range 7.5 M ≤ M1 ≤ 100 M and statistically weigh the systems by a Kroupa (2001) initial mass function. ZAMS mass ratios are selected at uniform intervals in the range 0.1 ≤ q ≤ 1 assuming a probability distribution f (q) ∝ q κ . In the fiducial simulation we assume a flat q distribution, i.e. κ = 0 (see e.g. Kuiper 1935;Sana et al. 2012;Kobulnicky et al. 2014). We draw orbital periods logarithmically-spaced in the range 1.41 days ≤ P ≤ 3.16 × 10 5 days with a probability distribution f (P ) ∝ P π . For M1 < 15 M we assume a flat distribution (π = 0) in log 10 P (Öpik 1924; Kobulnicky & Fryer 2007) while when M1 ≥ 15 M we assume a power law distribution in log 10 P with a slope π = −0.55 . The lower limit on log 10 P is also derived from , while we choose an upper limit sufficiently large to include non-interacting binaries as effectively single stars (de Mink & Belczynski 2015). The orbits are always assumed to be initially circular, since close orbits are expected to circularize due to tides or mass transfer before the first CC event (Belczyński & Bulik 1999;Hurley et al. 2002, although see also Eldridge 2009). For all stars in the fiducial run we assume a metallicity of Z = 0.02 (Anders & Grevesse 1989).
We calculate Roche lobe overflow (RLOF) mass transfer rates using the algorithm of Claeys et al. (2014). In this fiducial simulation, the mass transfer efficiency βRLOF -i.e. the relative ability of an accretor with mass Macc to accept mass lost by a donor at a rateṀ don -is defined as whereṀKH,acc=Macc/τKH, τKH is the Kelvin-Helmholtz timescale, and we set the parameter σ=10 (see Hurley et al. 2002;Claeys et al. 2014;Schneider et al. 2015). The matter that is not accreted during the mass transfer leaves the system. We assume it carries with it the specific angular momentum h of the accretor, i.e. h = γRLOFJ orb /(M don + Macc), where J orb is the total angular momentum, M don is the donor mass, and the loss parameter γRLOF = M don /Macc (see Soberman et al. 1997;van den Heuvel et al. 2017).
If the accretor star is sufficiently low-mass compared to the donor star, the system enters a phase of common envelope evolution (CE; Paczynski 1976). This critical mass ratio qcrit depends on the evolutionary stage of the donor star. We choose qcrit = 0.65, qcrit = 0.4 and qcrit = 0.25 when the donor star is a main sequence (de Mink et al. 2007), Hertzsprung gap, and red supergiant star respectively. We note that this last qcrit value might be relatively low (Claeys et al. 2014), i.e. the number of CE events with giant donors might be overestimated. However, R19 found no large variations in the number of runaway stars when increasing this value.
Once common envelope evolution begins, we model the evolution of the system using the αCEλ formalism (see e.g. Webbink 1984;de Kool 1990;De Marco et al. 2011), wherein a fraction αCE of the change in orbital energy ∆E orb from the inspiraling binaries is used to eject the envelope. For the binding energy parameter λ we use the analytic fit to the λg values of Dewi & Tauris (2000), which do not include the thermal energy of the envelope and its potential energy due to ionization. We assume in our fiducial simulation that αCE=1, i.e. all the liberated energy from the orbital inspiral is transferred to the envelope, which escapes to infinity at precisely the escape velocity from the system. Values αCE > 1 can be used to mimic the inclusion of other sources of energy (such as accretion luminosity, recombination energy, or nuclear burning) to aid the ejection of the envelope (e.g., Ivanova 2002;De Marco et al. 2011;Ivanova et al. 2013a).
We model the reaction of a system to a CC event following Tauris & Takens (1998). We assume that supernova ejecta leave the exploding star instantaneously since the ejecta speed is much larger than the orbital velocity of the binary. The mass loss from the companion due to stripping and ablation of its envelope as a result of the ejecta impact is treated following fits to the simulations of Liu et al. (2015) with a companion mass of M2 = 3.5 M .
The distribution of radio-pulsar proper motions (Shklovskii 1970;Gunn & Ostriker 1970) was the first piece of indirect evidence suggesting that asymmetries in supernova explosions could impart "natal kicks" to neutron stars at formation (Shklovskii 1970). Both asymmetric neutrino fluxes (Woosley 1987;Socrates et al. 2005) and large-scale density and velocity asymmetries in the star pre-collapse (e.g., Janka & Mueller 1994;Burrows & Hayes 1996) have been invoked to explain this natal kick, with the hydrodynamic-induced kick explanation currently more favoured over the asymmetric neutrino emission explanation (Holland-Ashford et al. 2017;Katsuda et al. 2018).
We draw kick magnitudes following Hobbs et al. (2005), who employ pulsar proper motions to infer that neutron star natal kicks follow a Maxwellian distribution with a root mean squared dispersion σ kick = 265 km s −1 . We draw natal kick directions isotropically in the frame of the collapsing star (see e.g. Wongwathanarat et al. 2013).
Since a large amount of ejecta can fall back onto the remnant following the CC, final natal kick velocities must be modulated accordingly, i.e. v kick → v kick (1 − f b ), where f b is the fallback fraction. We calculate fallback fractions depending on the core mass of the collapsing star following Eq. 16 of Fryer et al. (2012, hereafter F12), which also sets the mass of the remnant following a CC event. This algorithm generally assumes small (large) fallback fractions when the remnant is a neutron star (black hole) and therefore large (small) kick velocities.
In the event the natal kick launches the remnant directly towards the main sequence companion, we assume that the two bodies coalesce if the remnant enters the envelope of the companion (Leonard et al. 1994), i.e. if the post-CC periastron distance is less than the companion star's radius. We calculate the post-CC periastron distance following Eq. 57 of Tauris & Takens (1998), accounting for the change in momentum of the companion due to the impact of the expanding supernova shell and the stripping and ablation of the outer layers of the star companion. We assume in this calculation a shell velocity of 8200 km s −1 and an escape velocity from the companion of 800 km s −1 (Tauris & Takens 1998). We fit to Table 1 of Wheeler et al. (1975) to determine ablation and stripping mass fractions, assuming the companion is an n = 3 polytrope. In all, the impact of the supernova ejecta on the companion is a relatively small contribution to its total post-CC velocity, 10 km s −1 (Tauris & Takens 1998;Liu et al. 2015).
In total we draw 50 M1's x 50 q's x 100 P 's for a total of 250,000 binary systems. As we also draw 20 natal kick velocities and directions as well, we end up with a mock population of ∼5,000,000 evolved binaries.

Optimised Model Initial Conditions and Physical Assumptions
In addition to the fiducial model we also run a simulation with initial conditions and physics tuned to enhance the occurrence of high velocity ejections of companions from disrupted binaries. We discuss in Sec. 5.1 whether these tunings are consistent with contemporary observations and theory. If the size of the currently-known sample of HRS candidates cannot be matched by even this ad-hoc scenario, the BSS mechanism can be conclusively ruled out as a significant provenance for Milky Way HRSs. When a binary system is disrupted, the ejection velocity of the companion star is on the order of its pre-CC orbital velocity (Blaauw 1961;Eldridge et al. 2011;R19). The philosophy behind many of our changes to the initial conditions and physical parameters are therefore a) to increase the companion's pre-CC velocity as high as possible, and b) to increase the rate of binary disruption. Increasing the companion's orbital velocity implies producing harder binaries, which decreases the ease with which binaries are ionized. Indeed, tighter binaries may preferentially merge or remain bound and we therefore attempt to enhance the occurrence of physical configurations that specifically avoid those fates. Stable mass transfer can lead to orbital widening as the system tends towards equal star mass. We suppress the efficiency of the accretor to grow mass by setting βRLOF = 0. To encourage orbital hardening we also set γRLOF = 1 to increase the amount of angular momentum sapped from the system by ejected mass compared to our fiducial setup.
We want to encourage the occurrence of common envelope that leads to orbital tightening but simultaneously avoid binary mergers. We achieve this by choosing αCE = 10 -the maximum of the αCE range suggested by Hurley et al. (2002) -and also choosing qcrit = 1. A CE efficiency above might be achieved by tapping extra energy sources, e.g. in the thermal motions and ionization state of the envelope (De Marco et al. 2011;Ivanova et al. 2013a,b, see also Sec. 5.1.2).
Another physical ingredient that can in principle lead to tighter binaries is metallicity. At fixed stellar mass, stellar radii are smaller in metal-poor stars (see e.g. Burrows et al. 1993;Pols et al. 1998). Smaller stars enter the mass transfer stage later and experience less orbital widening, leading to faster orbital velocities (R19). In the optimized simulation we set the total stellar metallicity to Z = 0.0063, one-half dex lower than our fiducial value of Z = 0.02.
Finally, to increase the frequency of binary disruption, we boost the Maxwellian kick distribution dispersion to σ kick = 1000 km s −1 and set all fallback fractions to zero. We also set both the ZAMS mass ratio and initial log-period distribution slopes to π = κ = −1 to enhance the frequency of configurations where the companion is carrying most of the orbital velocity at ZAMS in a tight binary.

Other Model Variations
Our predictions depend on initial conditions for the binary systems and parametrized assumptions for binary evolution physics. In addition to the fiducial and optimized simulations, we run simulations where we vary the free parameters one-by-one. 2 . We set one parameter to its fiducial value while keeping all others to their optimized value to explore the impact each parameter has on the rate of high-velocity star ejections. To explore the possibility that the low-mass cores of less massive neutron stars experience smaller natal kicks, as suggested by some authors (e.g., Katz 1975;Arzoumanian et al. 2002;Pfahl et al. 2002;Podsiadlowski et al. 2004;Vigna-Gómez et al. 2018), we additionally perform a run with all parameters set to their optimized value except the natal kick is drawn from a double-peaked Maxwellian distribution similar to Vigna-Gómez et al. (2018) (see also R19). Following Pfahl et al. (2002) and Podsiadlowski et al. (2004), we draw from a kick distribution with σ kick,low = 30 km s −1 for remnants less massive than 1.35 M (Schwab et al. 2010;Knigge et al. 2011) while the kick distribution for more massive remnants remains at σ kick,high = 265 km s −1 . Table 1 summarizes the initial conditions and parameters assumed in each simulation.

Comparisons to Data
Each simulation contains a number N f of systems which eject a companion at a high velocity, here defined as vej ≥ 400 km s −1 . Properly accounting for the probabilities of each system and the finite stellar lifetimes of the ejected companions, we convert this to a number N f,now of stars in the Galaxy today with the same cut in ejection velocity. We estimate the number of systems -consisting of either isolated stars or binaries -that are formed per unit time in the Milky Way as SF R / M : we assume a constant star formation rate of SF R = 3.5 M yr −1 and the probability-weighted mean mass of a system is calculated as where f (M1) is a Kroupa (2001) initial mass function, f (q) ∝ q κ is the distribution of ZAMS mass ratios, and where we assume a binary fraction of 0 for M1 < 2 M and 1 for M1 ≥ 2 M . Eq. 2 gives M = 0.64 M for κ = −1 and 0.68 M for κ=0. Therefore, the birth rate of systems is SF R / M ≈ 5.1 − 5.5 yr −1 . A fraction F f of those systems eject fast stars and we calculate it by summing over the statistical weights pi of each sampled system that produces a fast star in our simulation given its initial configuration: These fast-ejected companions survive for a probabilityweighted average time t flight , where ∆t left,i is the remaining main sequence lifetime of the i-th fast-ejected companion. Putting this all together, the number N f,now of BSS-ejected fast stars in the Milky Way at any given time in the whole sky is therefore To roughly estimate a number N f,obs of fast companions currently observable in the whole sky, we cap ∆t left at 100 Myr under the assumption that fast stars with ∆t left > 100 Myr will be too far (and therefore too faint) to be detected at their current distances. This cap at 100 Myr is chosen to match more or less the maximum travel time seen in the current sample of hyper-runaway star candidates (see Sec. 3).
An alternative method to estimate N f,now would be to assume a CC supernova rate of ∼1.9 per century in the Milky Way (Diehl et al. 2006). Assuming all and only M > 8 M stars undergo CC supernovae, we can re-express F f as a fraction only of M1 > 8 M systems. Doing this and swapping SF R / M for the CC supernova rate in Eq. 5 yields very similar estimates for N f,now . Table 1. Summary of the input parameters to our binary synthesis models -the zero age main sequence (ZAMS) mass ratio distribution slope (κ), the ZAMS log-period distribution slope (π), the Roche lobe overflow (RLOF) mass transfer efficiency (β RLOF ), the RLOF angular momentum loss parameter (γ RLOF ), the critical mass ratio for common envelope evolution (q crit ), the common envelope efficiency (α CE ), the RMS dispersion of the natal kick Maxwellian distribution (σ kick ), the post-kick fallback fraction (f b ) and the total stellar metallicity (Z). See Sec. 2 for details.

SAMPLE OF MILKY WAY HYPER-RUNAWAY STARS
Our sample of Milky Way hyper-runaway candidates is derived from the Open Fast Stars Catalog 3 (OFSC; see Boubert et al. 2018), constructed using the AstroCats framework (Guillochon et al. 2017). The OFSC is a curated catalog of ∼500 hyper-velocity star candidates of all types, combining Gaia astrometry and photometry with supplementary properties such as stellar types and radial velocities from the literature. Spectra from SDSS (Abolfathi et al. 2018) or LAMOST ) are provided for each object when available. We extract full 6D astrometric measurements (position, proper motions, spectroscopic distance and radial velocity) from the OFSC for the 68 high-velocity star candidates identified as main sequence stars of spectral type A or B. For spectral types we refer only to the most recently published type for each star. For some hyper-runaway or hyper-velocity stars identified via colour-colour or T eff vs. log[g] selection cuts, it can be difficult to distinguish between early-type main sequence stars and lower-mass blue horizontal branch (BHB) stars since they populate similar regions of the Hertzsprung-Russel diagram (see e.g., Brown et al. 2006;Heber et al. 2008a;Brown et al. 2009). Projected rotational velocities and helium abundances derived via high-resolution spectra can break this degeneracy (see Heber et al. 2008a;Irrgang et al. 2018;Hattori et al. 2019). In cases where this ambiguity exists, we assume the star to be main sequence. Gaia proper motions are used for every star, though for HVS1, HVS10, HVS12 and HVS13 ) we additionally use the Hubble Space Telescope (HST) proper motion measurements of Brown et al. (2015), as they are more precise than the Gaia measurements even when ±0.5 mas yr −1 is added in quadrature to the uncertainties as suggested by Brown et al. (2018). HST proper motions are also used for the possibly-bound star B711 ), for whom the discrepancy between HST and Gaia proper motion measurements is large (see Irrgang et al. 2018). We remove HVS11 and HVS14 ) as well as HVS23 , for whom no proper motion measurements exist in the literature. When multiple distance or radial velocity measurements are provided, we take the value most recently published as of writing. For the 9 stars without explicit radial velocity uncertainties in Brown et al. (2009Brown et al. ( , 2012, we assume the reported ±11 km s −1 average uncertainty. We do the same for LMST HVS12 and LMST HVS19, assuming ±13km s −1 as the reported average uncertainty in Zhong et al. (2014). We keep sky positions for each star fixed and perform Monte Carlo (MC) resamplings over the remaining astrometric parameters, assuming Gaussian distributions and accounting for correlation between the proper motion components.
We integrate 1000 MC realizations for each of our 67 Aand B-type hyper-runaway candidates backwards in time up to a maximum of 1 Gyr at a fixed 0.1 Myr timestep in the MWPotential2014 potential outlined in Bovy (2015), a threecomponent potential model consisting of a spherical power law bulge potential, Miyamoto-Nagai disc (Miyamoto & Nagai 1975) and spherical NFW halo (Navarro et al. 1996).
3 https://faststars.space/ MWPotential2014 assumes a circular velocity at the Solar position of 220 km s −1 and a distance between the Sun and the Galactic Centre of 8 kpc. We assume local standard of rest UVW velocities of (11.1, 12.24, 7.25) km s −1 (Schönrich et al. 2010) and a height of the Sun above the stellar disc of 25 pc (Bland-Hawthorn & Gerhard 2016). We do not consider uncertainties in the Solar position or velocity, though we verify that these do not meaningfully affect our results. For the orbital traceback of each MC realization we determine the location of the last disc crossing, i.e. the (xGC, yGC) position when zGC = 0 in a Galactocentric Cartesian frame. For each realization we also track t f -the flight time from the disc crossing location to the observed position -and the ejection velocity vej in both the Galactocentric Cartesian frame and in the frame of the rotating disc. Ejection velocities in the corotating frame are computed using the rotational velocity curve associated to the MWPotential2014 potential. For MC realizations dynamically bound to the Milky Way, we record the location, velocity and flight time of each and every disc crossing up until the estimated main sequence lifetime of the star.
For each star in our sample we compute the probability P disc of being ejected from the Milky Way disc by taking the fraction of MC realizations which cross the disc at a Galactocentric distance 1 kpc ≤ rGC ≤ 25 kpc, imposing the cut at 1 kpc to remove stars ejected from the Galactic Centre and taking 25 kpc as the edge of the stellar disc (Xu et al. 2015). We also compute the probability P v>400km/s of a 'fast' ejection by taking the fraction of MC realizations which are ejected in the corotating frame at ≥400 km s −1 . This cut at 400 km s −1 is commonly cited as a limit for classical ejection mechanisms (see Irrgang et al. 2018, 2019, andreferences therein). This is a somewhat arbitrary velocity cut but our conclusions on the nature of the observed HRSs do not depend on the choice of the velocity threshold as long as it remains above a few hundreds km s −1 .
We select our early-type HRS candidates as those stars in our sample for which P disc >70 per cent and P v>400km/s >70 per cent. We additionally require that the Galactic Centre does not lie within the 1σ contour of the disc crossing location distribution. Our sample with these cuts applied results in 13 stars. The names, astrometry and spectral types of these stars can be seen in Table 2. We show as well the probability P ub of each HRS candidate to be unbound from the Galaxy, taken as the fraction of MC realizations whose total Galactocentric velocities exceed the MWPo-tential2014 escape velocity at their position. We also include in our sample as a special exception the oft-cited HRS candidate HD 271791, whose high space velocity, likely disc origin, and α-element enhancement are taken as evidence for a BSS ejection (Przybilla et al. 2008b). Note, however, that the natal kick velocity required in this case must be very large (Gvaramadze 2009). With an ejection velocity of 390 +70 −30 km s −1 , it just barely fails our P v>400km/s >70 per cent cut. Table 2. Observed properties of HRS candidates with a probability >70 per cent of being ejected from the Milky Way disc at a velocity > 400 km s −1 in the corotating frame. P ub gives the probability of each star being unbound from the Milky Way (see text). v GC is the star's total velocity in a Galactic Cartesian reference frame. Stars above the solid separating line only cross the Galactic disk within the nominal stellar lifetime of a star of their mass.

Name
Spec. Type  Table 3. Disc-crossing properties of 14 HRS candidates with a probability >70 per cent of being ejected from the Galactic disc at a velocity >400 km s −1 in the corotating frame. Uncertainties represent 1σ error intervals. Stars below the solid line cross the Galactic disc more than once within the estimated lifetime of a star with their mass. For such stars we only show the disc crossing properties for the most recent crossing. v ej,corotating denotes ejection velocities in the corotating frame of the Milky Way disc. P disc denotes the probability that the star is ejected from a Galactocentric distance 1 kpc ≤ r GC ≤ 25 kpc. P v>400km/s denotes the probability of that the star is ejected at v ej,corotating ≥ 400 km s −1 .  The Galactocentic Cartesian (yGC, yGC) locations of the most recent disc crossing for each MC realization of each star can be seen in Fig. 1. They range from stars ejected from a tightly-constrained region of the disc (e.g. B733) to stars with wide spreads of possible birthplaces (e.g. SDSS J115245.91-021116.2). In Fig. 2 we show the distributions of corotating ejection velocities for each hyperrunaway candidate. 10 of our 14 stars only cross the Galactic disc once within their estimated main sequence lifetime. At least 16 per cent of the MC realizations for two stars (SDSSJ081828.07+570922.1 and LMST HVS9) cross the Galactic disc twice within their MS lifetime, while the remaining two (LMST HVS19 and LMST HVS24) cross the disc more than twice. The possible ejection locations, velocities and flight times for the first crossing for each star are summarized in Table 3. With the exception of LMST HVS9 Luo et al. 2016) which is currently on an inbound orbit, every star in our sample has crossed the Milky Way disc within the last ∼100 Myr.
This sample of HRS candidates is neither volume-nor magnitude-limited. Since our sample is drawn from a variety of surveys and targeted searches with different sky coverage and completeness limits, we do not attempt to correct for these biases to determine a robust, corrected estimate for the number of HRS candidates in the Milky Way. Rather, this sample serves as a useful low-end estimate against which we compare the number of HRSs predicted by our binary population synthesis simulations.

RESULTS
We select from our simulations the main sequence companions ejected from disrupted binary systems, i.e. from systems which did not undergo a merger and did not stay bound after the CC of the primary. We further remove those systems which experience a collision between the remnant and the companion star post-CC (see Sec. 2.1). We show in Fig. 3   . Distribution of ejection velocities in the corotating frame from the Milky Way disc for 1000 Monte Carlo (MC) realizations of 14 HRS candidates. Different colours denote separate disc crossings, starting with the most recent. Shown for each star are P v>400km/s , the fraction of MC realizations where the star is ejected from the Milky Way disc during the most recent disc crossing at >400 km s −1 in the corotation frame, and P disc , the fraction of MC realizations where the star is ejected during the most recent disc crossing from a Galactocentric distance 1 kpc ≤ r GC ≤ 25 kpc. See Table 3 for more disc-crossing properties for each star.
both the fiducial (Sec. 2.1) and optimized (Sec. 2.2) binary evolution models, with the other models described in Sec. 2.3 not shown for brevity. We include for reference the ejection velocities (in the corotating disc frame) and stellar masses of the 14 stars in our sample of observed HRS candidates. Notice in the fiducial model (left panel) that the bulk of main sequence companions are ejected at ∼10 km s −1 regardless of mass, as reported in R19. There is in fact only a single main sequence companion in the fiducial model that is ejected at ≥ 400 km s −1 . By Eq. 5, this model predicts N f,now = 8 HRSs currently in the Milky Way and one BSSejected HRS currently observable, clearly at odds with the current known population of HRS candidates without even accounting for magnitude limits and sky coverage.
For the optimized model shown in the right panel of Fig. 3, notice that although the main peak of the vej distribution is still at ∼10 km s −1 , there is a second peak near vej ≈ 100 km s −1 with a high-velocity tail extending past our cut at vej = 400 km s −1 . The probability F f described by Eq. 3 of ejecting a fast companion is 2.5×10 −3 in the optimized simulation, ∼1.6×10 4 times more likely than in the fiducial simulation. This optimized simulation predicts N f,now 1.3×10 4 HRSs in the Milky Way, 3400 ejected less than 100 Myr ago -many more than the 14 HRS candidates in our observed sample. These typically have masses greater than ∼2 M , with a broad peak around ∼10 M that extends to the highest masses we probe. We point out that the BH former companions of these stars are even more massive, with a pre-CC mass distribution peaking at ∼30 M . The typical progenitor binaries of fast stars are therefore quite massive overall. We come back to this later in this section.
The results for the fiducial and optimized simulations as well as all other simulations described in Sec. 2.3 are summarized in Table 4. We show the probability F f of systems in the simulation to eject a companion at vej ≥ 400 km s −1 (see Eq. 3) along with the predicted number N f,now of vej ≥ 400 km s −1 stars in the Milky Way today and the number N f,obs with flight times less than 100 Myr (see Eq. 5).  Table 4. Summary of the frequency of fast ejections seen in the simulations described in Sec. 2 and summarized in Table. 1. We include the per-system probability F f of a ≥ 400 km s −1 ejection, the number N f,now of v ej ≥ 400 km s −1 hyper-runaway stars these simulations predict to be in the Milky Way today (see Eq. 5), the predicted number N f,obs of stars with a flight time less than 100 Myr, and the predicted fraction f RW 15 of all M > 15 M stars which are v ej > 30 km s −1 runaways (see text).

Model
F These results are useful in assessing which physical parameters and assumptions in our binary synthesis models most strongly affect the rate of high-velocity ejections. Keeping the notation of R19, we also include for reference the fraction f RW 15 of all M > 15 M stars predicted by each simulation to be runaway stars (vej ≥ 400 km s −1 ), properly accounting for the per-system probability and the duration of each evolutionary stage. This cut at 15 M corresponds roughly to O-type stars and makes f RW 15 a useful prediction against which to compare observations of Galactic runaways. In agreement with R19 we find that f RW 15 is robust against model variations and is limited to a few per cent or lower.
We defer a deeper discussion of O-type runaway fractions to Sec. 5.2.
From Table 4, it is clear that the frequency of fast star ejections is governed by only a few parameters. Setting the RLOF angular momentum loss parameter γ or the RLOF mass transfer efficiency β back to their fiducial values does not significantly change the probability for a companion to be ejected at a high velocity. The initial stellar properties, instead, are somewhat responsible for producing high velocity ejections. Notice that a reversion to our fiducial prescription for the ZAMS period power law slope (π = −0.55 for M1 > 15 M ) actually increases the probability to eject fast companions relative to the optimized simulation (π = −1 for all M1). This is perhaps counter-intuitive, as the steeper PZAMS slope in the optimized simulation should result in a larger probability for the low-PZAMS, high-v orb systems most likely to eject fast companions. However, it also significantly increases the probability of the PZAMS 10 day systems most likely to undergo a merger (see Sec. 4.1). Thus the slightly shallower PZAMS log-slope provided by the fiducial prescription ends up being more effective at encouraging the low-PZAMS HRS progenitor systems while not overproducing merging systems. Replacing our optimized PZAMS prescription with the fiducial prescription in the optimized simulation and all 'Opt. BUT' simulations does not change the qualitative results of this study.
Two other stellar initial condition parameters which affect the probability of a fast ejection are the ZAMS mass ratio distribution slope κ and total stellar metallicity Z. Increasing the total stellar metallicity back to its fiducial value of Z = 0.02 decreases N f,now by 50 per cent relative to the optimized simulation and returning κ from its optimized value of -1 to its fiducial value of 0 reduces the expected number of observable fast stars by 66 per cent.
The most impactful parameters controlling the prevalence of high-velocity ejections are the prescriptions for common envelope evolution and post-CC natal kicks. A phase of common envelope evolution appears indeed to be a common and crucial phase in the evolution of binaries that produce fast ejections: with the efficiency αCE set back to its fiducial value of unity, N f,now is reduced by 97 per cent. However, qcrit -also related to the CE evolution -does not significantly affect N f,now when set back to its fiducial value of 0.4. This is because systems capable of ejecting a main sequence companion at vej ≥ 400 km s −1 are already biased towards small mass ratios to maximize the orbital velocity of the companion. 94 per cent of systems ejecting a companion at vej ≥ 400 km s −1 in the optimized simulation have an initial q < 0.4.
With the kick dispersion σ kick set back to its fiducial value of 265 km s −1 (but fallback ratios still set to zero), the number N f,obs decreases by 99 per cent relative to the optimized simulation. Implementing the double-peaked Maxwellian natal kick distribution described in Sec. 2.3 decreases N f,obs by a similar amount. Mediating the natal kicks by returning to the F12 prescription for ejecta fallback fractions decreases N f,obs by 81 per cent. Together, these results demonstrate the relevance of black hole kicks for fast ejections of companions -black holes constitute 96 per cent of the remnants left behind by CC events in systems ejecting a companion at vej ≥ 400 km s −1 in the optimized simulation. The F12 prescription assumes large fallback fractions for black holes and thus very small natal kicks.
Small ZAMS mass ratios allow for larger companion orbital velocities and low stellar metallicities result in small stellar radii and therefore tighter orbits, but these contributions are minor when compared to the role of natal kicks and CE evolution. Included in Table 4 is a simulation with σ kick , αCE and f b reset to their fiducial prescriptions. Relative to the optimized simulation, the number N f,obs of t flight < 100 Myr stars observable in the sky is reduced by over 99.9 per cent and clearly demonstrates the importance of these parameters. The impact of natal kicks and CE evolution is also independent of κ and Z. We include in Table   4 a simulation with all of κ, π and Z set to their fiducial values. Relative to this simulation, additionally returning αCE, f b or σ kick to their fiducial prescriptions still results in a massive reduction in the number of companions ejected with vej ≥ 400 km s −1 . Note as well that the 'Opt. BUT αCE fiducial' and 'Opt. BUT σ kick fiducial' simulations predict a number of t flight <100 Myr hyper-runaway stars similar to our sample of 14 observed HRS candidates, and are thus not likely to be consistent with observations when accounting for completion and sky coverage limits. This indicates that all of αCE, σ kick and f b must be set significantly far from their fiducial prescriptions for the BSS to be consistent with current observations, as any single one set individually to its fiducial value is inconsistent.

The properties of progenitor binaries to fast ejections
The results presented above motivate the following picture for the production of an HRS through isolated binary evolution. Massive binaries with quite unequal mass ratios achieve small pre-CC orbital periods via CE evolution. As a result, the companion to obtains a very high pre-CC orbital velocities. RLOF evolution does not appear to play a role for the formation of HRS. The massive primaries undergo corecollapse events, leaving black hole remnants. Subsequently, strong natal black hole kicks are required to disrupt the binaries. The main sequence companion is then ejected with a velocity similar to its pre-CC orbital velocity. In Fig. 4 we illustrate the importance of tight binary orbits. In this plot we show how binaries in both the fiducial (left) and optimized (right) simulations populate the period pre-CC vs. the period at ZAMS plane. By the time of the first CC event, systems above the solid white line have had their orbits widened, generally via mass loss from stellar winds for large PZAMS and wind mass loss plus RLOF at shorter PZAMS. Systems below the solid line, on the other hand, have experienced orbital tightening, mainly via common envelope evolution. Systems at PpreCC = 0 in this figure have merged prior to the first core-collapse event. The contours in the right panel show the distribution only for systems which eject a companion at vej ≥ 400 km s −1 . These systems begin with periods of ∼10 1 -10 2 days at ZAMS and experience aggressive orbital shrinking due to common envelope evolution, reducing their orbital periods to ∼1 day. Conversely, the 10 PZAMS 100 day systems in the fiducial simulation (left panel) either experience orbital widening or have their orbits shrunk to the point of coalescence and experience a merger. Regardless, it is clear that the fiducial model does not possess the same reservoir of PpreCC∼1 day systems from which to eject companions at large velocities.
As the common envelope efficiency αCE is increased, less orbital tightening is required to fully eject the common envelope. It may seem counter-intuitive, then, to increase αCE in the optimized simulation, as it results in wider post-CE systems. In this study we are concerned with the very high tail of the velocity distribution of ejected companions, ejected from systems which may have otherwise merged if αCE were set lower. The priority for setting the common envelope efficiency in the optimized simulation is allowing systems to tighten while avoiding a merger. This is clear in Fig. 4, where systems at much higher PZAMS experience  . Distribution in the pre-core collapse period P preCC vs. zero age main sequence period P ZAMS plane for binary systems in the fiducial (left) and optimised (right) simulations. Bins are coloured by the total probability of all systems within the bin, scaled relative to the most likely bin. Solid white line shows a 1:1 correspondence. Systems lying above the white line experience orbital widening via wind mass loss or Roche lobe overflow mass transfer. Systems below the line experience orbital tightening through common envelope evolution. Systems at P preCC = 0 are systems that undergo a merger before the first core collapse event. White contours in the right panel show the 10%/30%/50%/70%/90% contours of the distribution of systems that eject companions at v ej ≥ 400 km s −1 a merger in the fiducial simulation than in the optimized simulation.
To further demonstrate the conditions and evolutionary stages required for an HRS ejection via the BSS mechanism, we present in Fig. 5 a schematic outlining the evolution of the most common HRS progenitor system in the optimized simulation. Weighting by the per-system probabilities and the remaining lifetimes of ejected hyper-runaway companions, at ZAMS (phase A) the typical binary consists of a massive (M1 ∼ 27 +14 −10 M ) primary and comparatively smaller (M2 ∼ 3.2 +2.2 −1.3 M ) companion on a shortperiod orbit (PZAMS ∼ 47 +40 −18 days). When the primary star first fills its Roche lobe (phase B), the companion cannot accept the material lost from the primary and the system enters a phase of common-envelope evolution. Dynamical friction within the common envelope tightens the binary, decreasing the orbital period to PZAMS ∼ 0.72 +0.51 −0.21 days. The envelope is eventually ejected and the core of the primary is exposed as either a naked helium-burning main sequence star or naked helium-burning Hertzsprung-gap star (phase C). At this time, the orbital velocity of the companion is already quite high (v orb ∼ 370 +50 −60 km s −1 ). When the primary undergoes a core-collapse supernova (phase D), the natal kick delivered to the remnant black hole is very large (v kick ∼ 1700 +600 −500 km s −1 ). This kick is sufficient to ionize the binary (phase E) -the companion is ejected with a velocity comparable to its pre-CC orbital velocity (vej ∼ 445 +87 −35 km s −1 ) while the BH is ejected at an even larger velocity, (vej ∼ 1500 +700 −500 km s −1 ).

DISCUSSION
We have shown in the previous section that the binary supernova scenario is unlikely to contribute significantly to the current population of observed HRS candidates, not unless the model prescriptions for a number of parameters are altered significantly from our fiducial prescriptions. In this section we determine whether these alterations are consistent with current understanding. To offer extra constraints on our models, we explore the predictions our fiducial and optimized models give for the number and kinematics of vej > 30 km s −1 runaway stars and of binaries that remain bound following the core-collapse of the primary. Finally we explore alternative mechanisms that may contribute to the existing population of Milky Way HRS candidates.

Binary evolution parameters that have the largest impact
We demonstrate in the previous section that the abundance of HRS ejected through the BSS scenario depends intimately on prescriptions for BH formation (natal kick and ejecta fallback fractions), common envelope evolution, and (to a lesser extent) the binary initial mass ratio distribution and stellar metallicity. In this subsection, we explore constraints on each of these prescriptions from the literature, assessing whether the altered prescriptions we assume in the optimized model are consistent -or at least not in conflict -with the current understanding of stellar binary evolution. We discuss the relevant parameters in 'chronological' order in the life of an evolving binary (see Fig. 5).

Initial conditions: the stellar metallicity and binary initial mass ratio
Of the relevant parameters we discuss here, a less impactful parameter which affects the frequency of vej ≥ 400 km s −1 companions is κ, the slope of the ZAMS mass ratio distribution. Though progress has been made in the past decade, attempts at constraining κ for high-mass binaries remain hamstrung by the difficulty of detecting even intermediatemass companions around massive stars (see e.g. Kobulnicky

A.
Initially tight binary Important parameters: metallicity, mass ratio B.  . Schematic depicting the evolution of a typical system in the optimized simulation which ejects an HRS. An unequalmass binary has a short period at zero age main sequence. Once the massive primary fills its Roche lobe, the system enters a common envelope phase, further tightening the orbit and increasing the velocity of the companion. When the primary explodes in a core-collapse supernova, the natal kick applied to the remnant black hole is strong enough to ionize the binary. The companion is ejected with a speed comparable to its pre-CC orbital velocity. & Fryer 2007;Duchêne & Kraus 2013). Our fiducial model adopts a flat distribution, i.e. κ = 0, consistent with observations of nearby Galactic O-type binaries  and O-type binaries in the Cygnus OB2 association (Kobulnicky et al. 2014). However, our choice of κ = −1 in the optimized simulation is also supported by several recent studies. Sana et al. (2013) investigate O-type binaries in the VLT-FLAMES Tarantula Survey of the 30 Doradus region in the LMC (Evans et al. 2011) and derive κ = −1.0 ± 0.4. Though they alter the functional form to allow for an excess of near-equal mass binaries, Moe & Stefano (2013) derive κ ≈ −1 in the Milky Way, SMC and LMC when studying eclipsing O-type and B-type binaries found in the Hipparcos (Lefèvre et al. 2009), OGLE-II (Wyrzykowski et al. 2004) and OGLE-III (Graczyk et al. 2011) data sets, respectively, though uncertainties remain large except in the case of the LMC. Even steeper mass ratio distributions have been determined by e.g. Dunstall et al. (2015), who investigate B-type binaries in the VLT-FLAMES Tarantula Survey and fit a mass ratio distribution with a slope κ = −2.8 ± 0.8. Therefore while our choice of κ = −1 in the optimized simulation is reasonable, larger and more complete samples are required to determine whether this accurately describes massive binary systems in the Milky Way.
The initial stellar metallicity also affects the frequency of HRS ejections. The effect of metallicity should not be overlooked, as evidence already exists of metal-poor stars ejected from the edge of the Milky Way disc (HD 271791;Heber et al. 2008b;Przybilla et al. 2008b). The total stellar metallicity in our fiducial simulation for all stars is the canonical Solar metallicity Z = 0.02 (Anders & Grevesse 1989) and one-half dex lower (Z = 0.0063) in our optimized simulation. While our fiducial Z is slightly higher than modern estimates (see Asplund et al. 2009, and references therein), the Sun does not appear to have an atypical metallicity when compared to other G dwarfs (Fuhrmann 2008;Holmberg et al. 2009) or B stars (Przybilla et al. 2008c) in the Solar Neighbourhood. Additionally, moderately-massive binary systems ejecting fast companions in the recent past would be expected to be more metal-rich than the Sun if anything, due to chemical enrichment over the past 4.5 Gyr. In light of this, assuming a metallicity of Z = 0.0063 over the entire Milky Way is difficult to justify. However, pockets of low metallicity in the Milky Way could contribute significantly to our current observed sample of HRS candidates. Additionally, the negative radial metallicity gradient of the Milky Way should not be neglected -although the stellar density drops off towards the edge of the disc, the stellar metallicity may reach more than 0.5 dex below Solar by the edge of the disc (see Lemasle et al. 2018, and references therein).

Late-stage evolution: the common envelope efficiency
More impactful than the ZAMS mass ratio or assumed stellar metallicity, an important physical parameter that has an governs the predicted number of HRSs is the common envelope efficiency αCE. R19 found that this has little impact on the bulk population and kinematics of runaways from the BSS considering all vej ≥ 30 km s −1 companions, since common envelope episodes are rather rare events in massive binary evolution (at least in comparison to stable RLOF mass transfer). However, here we focus on a rare outcome of massive binary evolution -the ejection of a companion at very large velocity. Therefore, a rare evolutionary channel can (and in this case does) play a large role in the formation scenario of the outcome of interest. To eject enough companions at large velocities to be consistent with current observations, we have shown that a common envelope efficiency significantly larger than unity must be prescribed, i.e. on top of the liberated orbital energy, other energy sources must be tapped to assist in ejecting the common envelope. These extra sources may include thermal energy and recombination within the envelope (e.g Han et al. 2002;Webbink 2008;Ivanova et al. 2013a), the enthalpy of the envelope (Ivanova & Chaichenets 2011) or nuclear fusion (Ivanova 2002;Podsiadlowski et al. 2010). Due to our presently-insufficient understanding of the physics involved in the common envelope phase, αCE is poorly-constrained and likely differs from system to system as the relevant timescales and energy sinks and sources vary (see Regős & Tout 1995;Zorotovic et al. 2010;Ivanova et al. 2013a). Estimates for αCE must be inferred from either observations of individual post-CE systems (e.g. Afs , ar & Ibanoǧlu 2008; Zorotovic et al. 2010), via binary population synthesis simulations (see Zuo & Li 2014, and references therein) or from detailed hydrodynamical simulations of individual systems (e.g. Sandquist et al. 1998;De Marco et al. 2011;Ohlmann et al. 2016;Fragos et al. 2019). While sample sizes remain low, several studies (e.g. Zorotovic et al. 2010;De Marco et al. 2011;Davis et al. 2012) identify systems where αCE > 1. The αCE 1 binary progenitors of HRSs however are rare systems, and observations of much more common systems may be only partially relevant. Whether αCE = 10 is possible in a rare and specific number of instances throughout the history of the Galaxy is currently observationally and theoretically unclear.

The end: CC natal kicks
Among all the physical parameters we explore, the post-CC natal kick and its mediation by ejecta fallback is most influential on the frequency of hyper-runaway ejections by the BSS scenario. We find in our optimized simulation that the primaries of the binary progenitors of HRSs nearly always leave a black hole remnant. This is an effect of HRS progenitor binaries maximizing the orbital velocity of the companion, trending towards large total masses and small mass ratios. Primaries in these progenitor binaries will almost certainly be massive enough to leave a BH remnant. Our fiducial treatment for the natal kick distribution for all remnants, however, is based not on BHs but on pulsar velocities. It follows Hobbs et al. (2005) who infer the distribution of 3D velocities for a sample of young (and therefore relatively unaffected by the Galactic potential) pulsars and fit a Maxwellian distribution with a root mean squared dispersion σ kick = 265 km s −1 (see also Lyne & Lorimer 1994). The strengths of black hole natal kicks are poorly-constrained when compared to pulsars (e.g. Repetto et al. 2012;Janka 2013;Mandel 2016;Janka 2017;Renzo et al. 2019b) -single BHs can only be detected via microlensing of background sources (Wyrzykowski et al. 2016) or accretion from the ISM (e.g. Fender et al. 2013). It is generally accepted that black hole natal kicks are similar to or lesser in magnitude than pulsar natal kicks (see Janka 2017;Atri et al. 2019;Chan et al. 2020). The σ kick = 265 km s −1 kick distribution for pulsars may already overestimate kick velocities, as several studies both theoretical and observational favour a double-peaked Maxwellian distribution with an additional low-velocity peak in the vicinity of σ kick = 30 km s −1 (e.g. Arzoumanian et al. 2002;Vigna-Gómez et al. 2018). In the fiducial model the natal kick is then mediated by the fallback of supernova ejecta as prescribed in F12 Eq. 16, i.e. v kick → v kick (1 − f b ), where f b is the fallback fraction. Our fiducial feedback prescription assumes large fallback fractions for BHs, therefore very small natal kicks effectively indistinguishable from zero (R19).
While a paucity of single BH observations continues to hamstring efforts to constrain BH natal kicks, surviving binaries containing a BH can offer some insight, at least on lower-velocity kicks insufficient to ionize the binaries. Evidence for non-zero BH kicks can be found in the Galactic latitude distribution and velocity distribution of black hole X-ray binaries (BHXRBs) (Repetto et al. 2012;Repetto & Nelemans 2015;Repetto et al. 2017), in particular the peculiar velocities of certain BHXRBs, e.g. GRO 1655-40 (Willems et al. 2005, XTE J1118+480 (Fragos et al. 2009) and GS 1354-64 (Atri et al. 2019. While at least some BHXRBs may have experienced kicks similar in magnitude to pulsars, recent works find BH kicks significantly larger than ∼ 100 km s −1 are not required to explain current BHXRB observations (e.g. Zuo 2015;Mandel 2016;Atri et al. 2019). Detailed 3D simulations of core-collapse supernovae accounting for the effects of ecjeta fallback find that kicks ∼ 100 km s −1 are possible in specific fallback scenarios (Chan et al. 2018;Chan et al. 2020) but favor more modest kicks in general. Naturally, systems which have experienced kicks strong enough to unbind the binary will of course not be visible as BHXRBs, therefore stronger kicks are not precluded by these studies. We explore in Sec. 5.2 whether comparing the velocities of BH-MS binaries in our simulations which remain bound post-CC to known BHXRBs can offer useful model constraints.
Additional constraints on BH natal kicks will improve with the determination of the mass distribution of runaway stars (R19) and further gravitational wave detections of binary black hole mergers (O'Shaughnessy et al. 2017;Wysocki et al. 2018). Non-zero BH kicks must be invoked to explain the misalignment between the orbital plane of the BH binary and the spin of the more massive companion (e.g. Kalogera 2000;O'Shaughnessy et al. 2017), which is encoded into the gravitational wave signal. On the other hand, the BH natal kick is bounded from above by binary BH merger rates and the mass distribution of runaway stars. Strong BH kicks would often disrupt massive binaries, skewing the mass distribution of runaways towards higher masses and sharply reducing coalescence rates for binary BHs from isolated binary evolution (e.g. Belczynski et al. 2002Belczynski et al. , 2016b. For a review of the formation and evolution of isolated compact object binaries, see Postnov & Yungelson (2014). The contribution of BH binaries formed not through isolated binary evolution but via dynamical assembly in dense stellar environments must be considered as well (see Benacquista & Downing 2013). On the basis of the observed LIGO binary BH merger rate, Wysocki et al. (2018) disfavor a BH kick distribution with σ kick 200 km s −1 (see also Belczynski et al. 2016a). To conclude, while the black hole natal kick prescriptions we assume in the optimized model cannot be directly ruled out due to a lack of observations of isolated single BHs, indirect evidence from pulsar kicks, BH XRB kinematics and BH-BH mergers generally do not offer support for this scenario.

Beyond HRS: Broader Implications of Model Variations
The binary evolution model variations we assume in this work impact not only the ejection rate of HRSs in general, but leave an imprint on other observables. In particular, the binaries in our simulations which remain bound post-CC can offer model constraints when compared to observations, as can the population of ejected stars with more moderate velocities (vej > 30 km s −1 ). In this subsection we remark on these populations.

On surviving binaries
We discussed in Sec. 5.1.3 how BH natal kicks can be constrained by BHXRBs and BH-BH mergers. Since we do not model X-ray emission via accretion and do not follow our systems beyond the first CC event, our systems cannot directly offer predictions for BHXRB and binary black hole populations. However, we can make inferences based on the abundance and kinematics of BH-MS binaries in our simulation which remain bound after the first CC. Even if they are not disrupted, the natal kick and the change in gravitational potential due to mass loss from the primary can impart peculiar velocities to these systems, possibly up to and exceeding runaway (vej > 30 km s −1 ) velocities (van Oijen 1989). Runaway stars still in orbit with a compact remnant could be identified as single-line (SB1) spectroscopic binaries. If the orbits are sufficiently tight to allow the black hole remnant to eventually accrete material from the companion, these systems may become visible as BHXRBs.
In the fiducial simulation, 8.7 per cent of systems remain bound as BH-MS binaries after the first CC event. Our simulations predict the systemic velocity vsys imparted to these systems, here defined as the velocity of the centre of mass of the post-CC BH-MS binary in the frame of reference of its pre-CC centre of mass (Kalogera 1996;Tauris & Takens 1998). In agreement with R19 (c.f. Fig. 11) we find that ∼80 per cent of the BH-MS systems in the fiducial simulation acquire vsys 0 km s −1 . This is due to our fiducial prescription for post-CC ejecta fallback following F12 -many BH progenitors experience complete ejecta fallback; the remnant experiences no kick at all and no mass is lost from the system. Only 4 per cent of the BH-MS binaries in our fiducial simulation achieve systemic velocities in excess of 30 km s −1 .
The optimized simulation, on the other hand, includes very different prescriptions for BH natal kicks and ejecta fallback, and therefore predicts a very different velocity distribution for surviving BH-MS binaries. Owing to the large disruption fraction, only 0.8 per cent of systems in the optimized simulation remain bound as a BH-MS binary after the first CC event. Of these, however, 99.9 per cent achieve vsys > 30 km s −1 and 95 per cent have a post-CC separation below 300 R . The optimized simulation therefore predicts runaway BHXRBs to be quite common. In fact, 45 per cent of BH-MS binaries achieve hyper-runaway speeds, vsys ≥ 400 km s −1 . The optimized simulation predicts that lone hyper-runaway stars outnumber hyper-runaway stars with BH companions by only a factor of 1.4. To date, no hyper-runaway BHXRBs or SB1 binaries have been observed. With a sample of 16 known BHXRBs, Atri et al.
(2019) find no BHXRBs that intersect the Galactic disc with a peculiar velocity above ∼ 200 km s −1 ; a systemic velocity achieved by 83 per cent of the BH-MS binaries in the optimized simulation. This absence of high-velocity BHXRBs or SB1 systems strongly disfavours our optimized model. It is worth mentioning, however, that without accounting for deceleration by the Galactic potential, BH-MS binaries in our optimized simulation travel for ∼17 kpc on average from their birth locations before the companion star leaves the main sequence. At such distances any X-ray emission may be difficult to detect and BHXRB samples would therefore be biased towards low-vsys systems remaining at low Galactic latitude.
Systems that remain bound after the first CC are unlikely to be disrupted by the eventual CC of the remaining MS companion (e.g., O'Shaughnessy et al. 2017) provided they do not merge beforehand. If the separation is sufficiently small, the resulting BH-BH binary may eventually become a gravitational wave source. Since we do not follow the evolution of our binary systems beyond the first CC event, our simulations cannot provide predictions for gravitational waves merger rates. However, from the different binary disruption fractions it is clear that the fiducial and optimized simulations would give different predictions for the rate and characteristics of BH-BH mergers. Therefore, upcoming statistical samples of BH-BH mergers will provide model constraints, albeit from a different part of the parameter space than what is relevant for the formation of hyper-runaway stars (see previous Section).

On massive runaway stars
Naturally, the model variations we explore to encourage the ejection of HRSs also strongly influence the fraction of systems which eject companions at 30 km s −1 < vej < 400 km s −1 and the velocity distribution of these 'runaway' stars. Properly accounting for initial configuration probabilities and the duration of each evolutionary phase, note in Table 4 that the fiducial simulation predicts that in total, 0.54 per cent of all O-type (M 15M ) stars in the Milky Way are runaways. The O-type runaway fraction increases to 1.6 per cent in the optimized simulation but never significantly exceeds ∼2 per cent in any simulation presented in Table 4. Although the likelihood of hyper-runaway ejections varies by four orders of magnitude among our chosen models, the O-type runaway fraction is much more robust against model variation. This consistent with R19 and the binary synthesis study of Eldridge et al. (2011) but in potential tension with observational findings that 10-20 per cent of O-type stars are runaways (e.g., Blaauw 1961;Tetzlaff et al. 2011;Maíz Apellániz et al. 2018) and the suggestion of Hoogerwerf et al. (2001) that ∼2/3 of runaway stars are due to disrupted binaries (though see Jilinski et al. (2010) for a challenge to this conclusion). Unless we have overlooked other mechanisms which greatly shrink binary orbits and eject companions at runaway speeds, the inability of even our optimized simulation to match the above observations suggests that the role of the dynamical ejection scenario has been under-estimated and/or observations of O-type stars are biased towards runaways.

Further caveats
There are caveats which should be considered when comparing the HRS populations predicted by the simulations to the current sample of observed HRS candidates.
Notice that we assume in Eq. 5 a constant Galactic star formation rate, a reasonable assumption for at least the last ∼2 Gyr (Snaith et al. 2014(Snaith et al. , 2015. Spatial or temporal variations in the global star formation rate in the Milky Way will affect the number of observed HRSs, especially the more massive and thus shorter-lived among these. A recent burst of star formation somewhere in the Milky Way disc will lead to an increase in runaway and hyper-runaway stars ejected from that location after a delay corresponding to at least the lifetime of the most massive stars we consider (∼50 Myr, Zapartas et al. 2017). Our assumption for the global Galactic star formation rate, 3.5 M yr −1 , is a common choice for binary population synthesis studies (e.g. Dominik et al. 2012), but rather optimistic when compared to recent results converging on 1-2 M yr −1 (Chomiuk & Povich 2011;Licquia & Newman 2015, and references therein). However, replacing the star formation rate seen in Eq. 5 with a lower value will reduce all N f,now values in Table 4 by the same factor and would not qualitatively change any results presented in this study regarding the conditions and physics required to eject fast stars via the BSS mechanism. A lower star formation rate in fact strengthens our conclusion that ejecting enough fast companions to be consistent with current observations of Milky Way HRS candidates requires prescriptions for black hole natal kicks, ejecta fallback and common envelope evolution that significantly differ from our fiducial assumptions for the overall binary population.
As mentioned in Sec. 5.1.2, one limitation of our approach is that we assume a single value of αCE applies to all binaries. In reality there is more likely a distribution of values, where even optimistically only for rare systems these are equal to our optimized value to αCE = 10. In the optimized simulation this same caveat applies to the total stellar metallicity Z and fallback fraction f b -we apply a single value to the entire binary population rather than allowing for multiple populations to be present as minority fractions of the broader population. A model where Z and αCE in particular are allowed to vary is a venue worth exploring in the future.
Finally, it is worth pointing out that the parameters chosen for our optimized simulation were chosen based on simple assumptions about their role in enhancing the number of short-period, high-orbital velocity systems. We do not fully explore the highly-dimensional space of initial conditions and binary evolution parameters to find the truly optimal set of prescriptions. Some of our assumptions may be overly simplistic, e.g. see Sec 4 where the fiducial PZAMS prescription is in fact more effective than the "optimized" prescription at ejecting fast companions. A more complete exploration of the entire parameter space and the correlations between parameters is a possible future direction for this study.

Alternative origins of hyper-runaway stars
Here we only consider main-sequence stars ejected from a disrupted binary system. If the BSS mechanism is unable to contribute significantly to the known population of Milky Way main sequence hyper-runaway candidates, which mechanism is responsible? One possible mechanism includes the dynamical ejection scenario (DES) wherein stars are ejected via three-or four-body interactions in dense clusters (Poveda et al. 1967;Leonard & Duncan 1990;Leonard 1991;Fujii & Portegies Zwart 2011;Banerjee et al. 2012). Oh & Kroupa (2016) explore the ejection of massive stars from moderately massive (M cl ≈ 3000 M ) star clusters under diverse initial conditions. Ejections in excess of 200 km s −1 are very rare except in the model where the cluster consists of only initially massive (M1 ≥ 5 M ) binaries on short, circular orbits. In a suite of similar simulations, Perets &Šubr (2012) infer a dynamical ejection rate for vej > 450 km s −1 B-type HRSs from young clusters of ∼10 Gyr −1 in their most realistic model, and conclude that the DES is unlikely to contribute significantly to the Galactic population of hyper-runaway or hyper-velocity stars. Since the majority of ejections from dynamical encounters happen in the very early evolution of a cluster (e.g. Oh & Kroupa 2016), observations of the stellar dynamics in very young regions (such as the Orion Nebula Cluster, e.g., Schoettler et al. 2019Schoettler et al. , 2020 in the LMC, e.g., Lennon et al. 2018;Renzo et al. 2019a) can constrain the DES scenario, including for the rare very fast ejections.
Beyond the BSS and the classical dynamical ejection scenario, possible but less-studied scenarios have been suggested. A first potential pathway is ejections which involve three-or four-body interactions involving massive stars. Gvaramadze (2009) point out three possible ejection channels: (i) the breakup of a an unstable hierarchical triple composed of massive stars (inner binary of two ∼50 M stars with a ∼10 M outer companion, (ii) binary-binary dynamical interactions where massive stars (∼20-40 M ) compose at least one member of the binaries, and iii) exchange encounters between hard, massive binaries and 200 M very massive stars, which form via runaway collisions of ordinary stars in the cores of young massive star clusters (Portegies Zwart & McMillan 2002). All three mechanisms are capable of ejecting massive companions far in excess of 400 km s −1 .
Second, encounters between stellar binaries and an intermediate-mass black hole (IMBH) in analogy to the Hills mechanism may also contribute to the population of highvelocity stars ejected from the Galactic disc or halo (Pfahl 2005;Gualandris & Portegies Zwart 2007;Hopman 2009;Sesana et al. 2012). Fragione & Gualandris (2019) assign realistic IMBH masses to the centres of Milky Way globular clusters and study the ejection of high-velocity stars. They infer a fast star ejection rate from this IMBH channel of ∼10 −4 yr −1 with an ejection velocity distribution which peaks at vej ≈ 300 km s −1 with a tail extending to ∼2000 km s −1 .
Finally, we must also consider the possibility that although the backward trajectories of observed HRSs point towards the Galactic disc, they may have been born elsewhere. Main sequence HRSs originally of extragalactic origins may serendipitously intersect the Milky Way disc and appear as disc-ejected HRSs. These stars could come from the tidal disruption of infalling dwarf galaxies as suggested by (Abadi et al. 2009) (Sherwin et al. 2008), though none of our HRS candidates seem to point from the LMC or M31. Any or all of the above origins and mechanisms could contribute to the Galactic population of hyper-runaway.

SUMMARY AND CONCLUSION
We have performed an extensive suite of numerical simulations to ascertain the probability with which disruptions of binary systems by core-collapse supernovae are able to eject main sequence companions at large velocities (vej ≥ 400 km s −1 ). For each simulation we convert this probability to a prediction for the number of these so-called 'hyperrunaway stars' (HRSs) in the Milky Way which should still be reasonably nearby with t flight < 100 Myr. We systematically vary parameters which dictate the initial conditions and physics governing the evolution of binary systems. In doing so, our aim is to rigorously test under a variety of models whether the binary supernova scenario contributes significantly to the observed population of early-type HRS candidates seemingly ejected from the Milky Way disc at high velocities. By varying parameters one-by-one, we determine which are most important when considering the ejection of companions at large velocities. Our findings can be summarized as follows: • In our fiducial model, which generates and evolves binary systems according to free parameters chosen to be consistent with contemporary observations and theory, corecollapse supernovae in binaries do not eject companions at large velocities often enough for the binary supernova scenario to contribute significantly to the known population of HRS candidates (see Fig 3).
• By varying relevant parameters one-by-one, we find that the probability of fast companions ejections depends most intimately on the common-envelope efficiency αCE and the parameters which dictate the magnitude of the post-CC natal kick imparted to remnant black holes, i.e. the ejecta fallback fraction f b and natal kick distribution dispersion σ kick (see Table 4). The αCE parameter indirectly governs the amount of orbital tightening -and therefore the increase in orbital velocity -that occurs during the common-envelope phase. A large orbital velocity in turn results in a greater ejection velocity for the companion if the system is disrupted. The natal kick parameters f b and σ kick , in turn, dictate the frequency with which binary systems are disrupted immediately following the core-collapse of the primary.
• To a lesser degree, the probability of fast companion ejection increases with decreasing total stellar metallicity and with an increasing proportion of highly-unequal initial mass ratios in the binary population. At fixed stellar mass, stars of lower stellar metallicity are smaller in radius and can therefore orbit at a smaller separation without undergoing a merger. At fixed system mass, a smaller mass ratio between the companion and primary results in a greater orbital velocity for the companion.
• While the rate of HRS ejections via the binary supernova scenario depends strongly on the above parameters, the runaway (vej > 30 km s −1 ) fraction among O-type stars is robust against model variations and is limited to 2 per cent.
• The prototypical progenitor binary of an HRS ejected by the binary supernova scenario is composed of a massive, unequal-mass binary on a short initial period. Upon the primary filling its Roche lobe, the system enters a phase of common envelope evolution and hardens further. The core collapse of the primary leaves a black hole remnant, which experiences an asymmetric post-CC natal kick extreme enough to disrupt the binary. The main sequence companion escapes with an ejection velocity similar in magnitude to its pre-core collapse orbital velocity (vej 400 km s −1 , see Figs. 3 and 5). The remnant black hole is ejected at a much larger velocity, vej 1700 km s −1 .
• For a model to predict a number of t flight <100 Myr HRSs comparable to the current observed population of HRS candidates, it is not sufficient to tune one of the aforementioned parameters in the direction which favours the ejection of fast companions. All of the common envelope efficiency αCE and the black hole natal kick parameters f b and σ kick must be tuned in favourable directions for our population synthesis model to predict a number of hyper-runaway stars consistent with observations (see Table 4).
• To match observations, the relevant parameters above must be tuned to values or prescriptions which are improbable or disfavoured -but not ruled out -by contemporary constraints on these parameters from the literature.
From the results summarized above, we may conclude that the binary supernova scenario is unlikely to contribute significantly to the current known population of HRS candidates. This conclusion is broadly consistent with other works exploring the binary supernova scenario (e.g. Portegies Zwart 2000; Eldridge et al. 2011;Renzo et al. 2019b). With our detailed quantitative investigation we outline for the first time the evolutionary path of hyper-runaway stars ejected via the binary supernova scenario and the characteristics of their binary progenitors. Future dedicated theoretical and observational works may be used to test the existence of such progenitors.