© P. Neunteufel et al. 2022

1. Introduction

For the purposes of measuring cosmological distances, supernovae (SN) of Type Ia, in their role as standard candles, occupy a dominating position in the observational toolset (leading to a Nobel Prize in 2011, see Perlmutter et al. 1999). The usefulness of transients of this type relies heavily on the well-established relationship between the peak-luminosity light curve evolution (Phillips 1993) of the event, allowing the observer to independently infer the brightness of the event from parameters accessible to observation. The observational characterisation, which includes spectral analysis (Blondin et al. 2012; Hillebrandt et al. 2013) of a SN Ia or other type has to be distinguished from theoretical classification of the production mechanism. For example, in the context of the presumed production mechanism, a collapse of the iron core of a massive star (core collapse SNE, CCSNe) may result in a variety of spectral classifications depending on the state of the progenitor in terms of, for example, its mass, its chemical composition, its structure and other parameters (see e.g. Janka 2012; Janka et al. 2016). In the same vein, SNe Ia may result from a variety of different progenitor scenarios. It is generally accepted that SNe Ia result from the thermonuclear explosion of a white dwarf star (WD) – the degenerate, burnt-out remnant of a relatively low mass (≲8 M_⊙) main sequence (MS) star (e.g. Weidemann 2000; Cummings et al. 2018). WDs are inert in a thermonuclear sense and, if left to their own devices, will simply cool radiatively until they reach thermal equilibrium with the surrounding medium. It is therefore further generally accepted that WDs require material to be donated from another source in order to reinitiate nuclear burning, resulting in an explosion. The process by which this reignition occurs and then proceeds through the WD is currently not sufficiently well understood (see e.g. Hillebrandt & Niemeyer 2000; Hillebrandt et al. 2013; Livio & Mazzali 2018; Ruiter 2020). While a number of different progenitor scenarios for SNe Ia have been proposed in recent decades, as will be elaborated on below, a clear favourite for the dominant progenitor channel has not been forthcoming.

The likely nature of SNe Ia as accretion induced thermonuclear explosions associated with WDs puts them into a conceptual family of thermonuclear transients which also includes classical novae (see e.g. José et al. 2006; Chomiuk et al. 2021), cataclysmic variables (CVs, see e.g. Gänsicke 2000) and AM Canum Venaticorum (AM CVn, which have been associated with SN Ia progenitor candidate systems, Podsiadlowski et al. 2003; Nelemans 2005; Solheim & Yungelson 2005; Solheim 2010; Piersanti et al. 2015) and its analogues. Additionally, a number of other transients originating from WDs have been either proposed or classified. Among these are SNe of Type Iax (Li et al. 2003; Foley et al. 2013) and other peculiar SN Ia variants (Woosley & Kasen 2011). All of these together, such as thermonuclear events, including ordinary SNe Ia, originating from WDs which are brighter and more energetic than CVs and classical novae and for the purposes of this paper, form the family of ‘thermonuclear SNe’. From this point onwards, we refer to the spectral classification of a SN of Type Ia only and exclusively as ‘ordinary’ SNe Ia. Attempts to identify the progenitor systems of observed individual ordinary SNe Ia have generally yielded ambiguous results. Counter-intuitively, investigations in this direction have actually been more successful for SNe Iax than for ordinary SNe Ia, McCully et al. (2014) having identified a luminous blue source at the location of SN 2012Z, commensurate with either a massive main sequence star or a hydrogen depleted star of a mass ∼1.0 M_⊙.

Progenitor scenarios for thermonuclear SNe are generally accepted to fall into two distinct categories, these being the single degenerate and the double degenerate channels. In the double degenerate channel, the progenitor binary is assumed to consist of two WDs (which may be of different mass and composition). Under the influence of angular momentum loss due to gravitational wave radiation (GWR Landau & Lifshitz 1975), and assuming a not exceedingly large initial orbital separation, this progenitor system will become tight enough within Hubble-time for the less massive WD to fill its Roche-lobe. Then, given a mass ratio not too different from unity, it will undergo dynamically unstable Roche-lobe overflow (RLOF). Unstable ignition of either the accreting WD, both WDs individually, or an object resulting from the merger of the two components, is expected, leading to a SN (see e.g Webbink 1984; Tornambe & Matteucci 1986; Cameron & Iben 1986; Sim et al. 2010; Pakmor et al. 2012, 2013). The number of candidate progenitors in this channel is surprisingly limited, with only two currently known candidates (KPD1930+2752, HD 265435, Maxted et al. 2000; Pelisoli et al. 2021 respectively). In the single degenerate channel, the accreting (degenerate) WD is assumed to receive mass from a non-degenerate companion of some description. The donor star, according to different sub-channels, may be an ordinary main sequence star, an evolved post-main sequence star, either of which would chiefly donate hydrogen-rich material, or a hydrogen depleted star which would chiefly donate helium enriched material (Nomoto 1980, 1982a,b; Wang et al. 2017). Again, a number of different ignition and explosion mechanisms have been proposed, a thorough discussion of which would extend beyond the limits of this article (though see Ruiter 2020 for an extensive review). Between these proposed mechanisms and a lack of direct observational evidence, the structure, in terms of component masses and stellar types involved, of the candidate progenitor binaries, as well as the identity of the spectral type resulting from a particular channel and sub-channel, remain obscured.

One feature distinguishing the single degenerate channel from the double degenerate channel is that the donor star is generally expected to survive the destruction of the accretor (Podsiadlowski 2003; Liu et al. 2013; Bauer et al. 2019; Neunteufel 2020; Liu & Zeng 2021; Liu et al. 2021). Therefore, in this channel, as the accreting WD is disrupted in a SN, the former donor star is flung away, under conservation of its pre-SN orbital angular momentum, from the centre of mass of the progenitor binary (see e.g. Hansen 2003; Justham et al. 2009). The velocity of this putative runaway star depends, as will be discussed below, on the star’s pre-ejection mass and structure, with highest runaway velocities in the single degenerate channel associated with the helium donor sub-channel.

In many cases, the ejected helium star will be fast enough to be unbound from the Galactic potential, in the process becoming a hypervelocity star (HVS). HVSs, of which the type under consideration here is not the only variant, have attracted considerable interest over the years. Specifically, higher mass and main sequence candidates are expected to form during interaction of a binary star with a supermassive black hole (Hills 1988; Brown 2015; Hamers & Perets 2017). Notably, Koposov et al. (2020) reported on the discovery of an object originating at the Galactic centre. HVS of Galactic centre origin have been utilised in probing the Galactic matter distribution (Yu & Madau 2007) and the environment of the Galactic centre (Brown et al. 2018; Generozov & Madigan 2020; Evans et al. 2022). This present work should be considered in the same line as others focused on the kinematic and population properties of such objects (e.g. Brown et al. 2009; Lu et al. 2010; Rossi et al. 2014), though, unlike these, considering ejected He-stars of SN origin (Geier et al. 2015a; Vennes et al. 2017; Shen et al. 2018). We mention that origins as former companion stars in core collapse SNe have also been considered (Renzo et al. 2019; Evans et al. 2020). In any case, HVS stars of all masses and configurations are a topic of active research (Raddi et al. 2021; Li et al. 2021) and observational campaigns focused on their detection, especially in the era of the Gaia Instrument (Gaia Collaboration 2016), are ongoing (Irrgang et al. 2018; Bromley et al. 2018; Raddi et al. 2021; Marchetti 2021; Irrgang et al. 2021) and not limited to the Galaxy (Evans et al. 2021). Most recently, ejection of high velocity stars from dynamical disruption of triple systems has also been considered (Hamers et al. 2022). We note that HVS are, conceptually, a variant of runaway stars – stars with a velocity excess over their direct association – as first considered by Blaauw (1961). In this study we focus on the single degenerate channel, with a particular focus on the helium donor sub-channel (hereafter: SD-He channel). Specifically, we attempt to resolve the following: firstly, the structure and distribution of a hypothetical population of runaway stars originating from the SD-He channel; secondly, how this structure depends on initial assumptions on the parameter space, specifically with respect to terminal accretor mass; thirdly, whether there are any candidate members of this population; fourthly, which conclusions on the structure of the progenitor binary can be drawn from the observational record; lastly, the theoretically predicted occurrence rates current observational record supports for this channel.

This paper is organised as follows: In Sect. 2 we expound on the context of this study and motivate our work. In Sect. 3, we discuss our physical assumptions and explain our numerical methodology. In Sect. 4 we summarise the kinematic properties of the population. In Sect. 5 we use the results presented in the previous section to qualitatively define criteria for adherence to this population and discuss possible sub-divisions. In Sect. 6 we discuss the sensitivity of population parameters to assumed initial conditions and define corresponding observables, broadly taking into account constraints on the observability of population members. In Sect. 7 we attempt to identify possible observed members of our predicted population and categorise them in their expected subpopulations. We comment on astrometric parameters of our population in Sect. 8. Using the identified members as well as the predictions resulting from our simulations, we estimate the event rate for SD-He ejections of US 708-like stars as supported by the currently available observational record as well as calculate the number of detections necessary to support these ejections as significant contributors to the Galactic SN rate in Sect. 9. We conclude and summarise our study in Sect.10.

2. Motivation and context

2.1. Thermonuclear explosions of WD stars put into context

As stated above, this article is chiefly concerned with survivors of the single-degenerate helium donor channel for thermonuclear SNe. In recent years, three major channels and a number of minor ones (in the sense of community interest) have emerged: firstly, the double degenerate channel, which has already been commented upon, which involves a progenitor system composed of two WD stars being driven to RLOF by emission of GWR. This channel is usually not supposed to produce a surviving remnant (Webbink 1984), though exceptions leading to the ejection of a fast runaway have been proposed (see discussion in Sect. 2.3); secondly, the hydrogen donor channel (also: classical mechanism) where the accreting WD receives mainly hydrogen-rich material from a companion star in on the main sequence or in a red giant configuration. In this scenario, the accreting WD is usually assumed to reach the Chandrasekhar mass, leading to carbon ignition at its core (see reviews by Hillebrandt & Niemeyer 2000; Wang & Han 2012; Ruiter 2020). Due to the size of the assumed donor star, this mechanism is incapable of producing a high-velocity runaway. Thirdly, the donor star is assumed to be a hydrogen-depleted helium star. The transferred material will thus be helium enriched. As shown by Kato & Hachisu (2004), Piersanti et al. (2014) and others, depending on the mass transfer rate (and assumed input physics, such as rotation), the transferred material will be accumulated quiescently (Ṁ ≲ 3 × 10⁻⁸ M_⊙ yr⁻¹), undergo violent flashes (3 × 10⁻⁸ M_⊙ yr⁻¹ ≲ Ṁ ≲ 1 × 10⁻⁷ M_⊙ yr⁻¹), weak flashes (1 × 10⁻⁷ M_⊙ yr⁻¹ ≲ Ṁ ≲ 5 × 10⁻⁶ M_⊙ yr⁻¹), or steady burning (Ṁ ≳ 4 × 10⁻⁶ M_⊙ yr⁻¹). If sufficient material is accumulated, the WD’s mass may grow to the Chandrasekhar mass, eventually undergoing detonation via the classical mechanism or, if helium ignition conditions are reached, produce a He-detonation, igniting the CO core at the He-CO boundary. The latter may happen at terminal accretor masses different from the Chandrasekhar mass.

2.2. Chandrasekhar mass explosions

For the purposes of this study, we need to distinguish the classical mechanism, and Chandrasekhar mass double detonations. The classical mechanism assumes that a WD accretes sufficient (hydrogen-rich) material from a companion to approach the Chandrasekhar mass. Close to the Chandrasekhar mass, the degenerate electron gas will approach the relativistic limit, leading to a change in the equation of state. The consequent increase in density precipitates thermonuclear ignition and explosion of the WD (Hoyle & Fowler 1960; Arnett 1969; Hansen & Wheeler 1969). Detailed modelling of the explosion has seen significant advances over the years (see e.g. Nomoto et al. 1984; Hillebrandt & Niemeyer 2000), with a number of different mechanisms having been proposed, but none having been identified conclusively as dominant. Further questions persist as to the ability of a system of this structure (i.e. a WD + hydrogen-rich companion) to provide the necessary conditions for the WD to grow to the Chandrasekhar mass and their contribution to the total SN Ia rate has been estimated to be < 5% by Gilfanov & Bogdán (2010), but see for example Di Stefano (2010). However, since, for hydrogen-rich MS stars, expected ejection velocities lie well below the velocity required to unbind the ejected star from the Galactic potential, the classical mechanism is not considered in this study.

Instead, for our purposes, we consider Chandrasekhar mass explosions to result from the same mechanism as sub-Chandrasekhar mass explosion (as discussed in greater detail in Sect. 2.3) and in the same type of system, such as He-star+WD binaries, likewise producing a runaway He-star. The possibility of a near-Chandrasekhar mass double detonation is usually predicated on inclusion of microphysics, such as magnetic torques, as proposed by Spruit (2002) and modelled by Neunteufel et al. (2017, 2019). However, Chandrasekhar mass explosions originating from He-accretion have also been connected to SNe of Type Iax and have been predicted to leave a partially burnt remnant (Gamezo et al. 2005; Branch et al. 2004; Phillips et al. 2007; Foley et al. 2009; Kromer et al. 2013, 2015; Fink et al. 2014; Magee et al. 2016; Lach et al. 2022).

2.3. Sub-Chandrasekhar mass explosions

Faced with difficulty in explaining how the future SN progenitor (i.e. the accreting WD) can accumulate sufficient mass to reach the Chandrasekhar mass, the question arises whether the star is actually required to reach it, leading to any number of sub-Chandrasekhar mass scenarios. The most well know of these was originally proposed in the 1980s by Nomoto (1980, 1982b), but has consistently generated interest in the intervening years (e.g. Livne 1990; Livne & Arnett 1995; Yoon & Langer 2004a,b; Fink et al. 2007; Sim et al. 2010; Woosley & Kasen 2011; García-Senz et al. 2018; Neunteufel et al. 2016, 2017, 2019; Gronow et al. 2021a,b; Shen et al. 2021). This mechanism, known as double detonation mechanism achieves detonation of the WD’s CO core by treating an accumulated, quiescent helium layer, donated by a hydrogen depleted companion, as a thermonuclear fuse. Here, compressive heating, additional nuclear heating via the NCO-effect (Hashimoto et al. 1986) and (depending on included microphysics, viscous heating), precipitate ignition and, under certain circumstances, detonation of the accumulated helium layer (the first detonation). The first detonation then propagates into the CO-core, leading to explosive carbon ignition (the second detonation). The amount of helium needed for the first detonation to occur is being debated in literature but is generally thought, depending on included microphysics and the temperature and mass of the accreting body, to lie in the range of 0.05 M_⊙ − 0.4 M_⊙. The minimum mass of the WD (prior to the terminal accretion episode) has been argued to be as low as 0.75 M_⊙ (Livne & Arnett 1995). One peculiarity of progenitors in this channel is the requirement for the gravitational wave merger timescale to be short enough for the He-donor, most likely a subdwarf B or subdwarf O star (Heber 2009, 2016; Preece et al. 2018, 2019), to fill its Roche lobe while still in its core-He burning phase (which is expected to last for ∼100 Myr). This requires the binary at formation of the terminal configuration of the system to be compact. Compactness is customarily argued to be a result of the formation of either the subdwarf BO star or the WD in a common envelope (Han & Podsiadlowski 2004; Han et al. 2002, 2003; Justham et al. 2011; Bauer & Kupfer 2021; Kruckow et al. 2021). Only two candidate progenitor systems in this channel have been characterised so far. One is CD-30°11223 (Vennes et al. 2012) and, potentially and very recently, PTF1J2238+7430 (Kupfer et al. 2022).

Another proposed mechanism in this category bears mentioning: The so-called prompt double detonation or dynamically driven double degenerate double detonation mechanism. Here, turbulent ignition of a thin < 0.05 M_⊙ helium layer, accreted during the dynamical disruption of a degenerate companion, for example a He-WD or hybrid HeCO WD (e.g. Pakmor et al. 2021), ignites a secondary detonation of the accretor’s CO core (Pakmor et al. 2013). This mechanism is especially important to keep in mind in the context of this study, as one of the candidates (LP 398-9, see Sect. 7) has been claimed as a potential survivor of this scenario, as argued by Shen et al. (2018), but see Meng & Podsiadlowski (2017).

2.4. Super-Chandrasekhar mass explosions

A straightforward scenario for the formation of a super-Chandrasekhar mass configuration that may lead to explosion is the merger of two sub-Chandrasekhar mass WDs with a combined mass larger than the Chandrasekhar mass. In a single-degenerate scenario, explosions well above the canonical Chandrasekhar mass may be achieved if the WD is stabilised against ignition or collapse by rotation. Models of rotating WDs stable far above the canonical Chandrasekhar mass were provided by, for example, Ostriker & Bodenheimer (1968) and Hachisu (1986). A distribution of the SN Ia explosion mass over a range of 0.08 M_⊙, as a consequence of solid-body rotation in the accreted material, was investigated by Uenishi et al. (2003) as a potential source of diversity in observed SNe Ia (see also Yoon & Langer 2004b). However Yoon & Langer (2005) found that, for plausible assumptions, the angular momentum gained by accreting WDs might allow differentially-rotating WDs which reach 2.2 M_⊙. When accretion is no longer able to maintain the necessary differential rotation, the WD may lose or redistribute spin angular momentum (e.g. under the influence of the r-mode instability or tidal interaction, see e.g. Yoon & Langer 2004b), and so ignite and explode. In the context of this paper, explosions of such highly super-Chandrasekhar WDs from single-degenerate SN Ia progenitors could increase the maximum ejection velocity of hypervelocity stars. However, depending again on the microphysics of the WD (see e.g. Neunteufel et al. 2017), it is debated whether, and under what circumstances, differential rotation can be maintained in accreting WDs (see, e.g. Yoon & Langer 2004b; Piro 2008; Hachisu et al. 2012). Nonetheless, even if differential rotation cannot be maintained in the accreting WD then, if ignition of the WD due to other factors (such as viscous heating) can be avoided, solid-body rotation may allow the WDs in some single-degenerate SN Ia progenitors to reach ≈1.5 M_⊙ without explosion or instability (e.g. Hachisu 1986; Uenishi et al. 2003).

Justham (2011) and Di Stefano et al. (2011) argued that considering the increased WD mass at explosion as a consequence of accreted angular momentum plausibly also resolves some of the arguments which had been used against single-degenerate SN Ia progenitors, with the latter paper describing the scenario as the “spin-up-spin-down” model (see also, e.g. Hachisu et al. 2012 for a clear investigation of potential advantages of this model). A common simple picture of this scenario is of the subset of cases in which mass accretion totally ceases before the explosion, with the WD then spinning-down towards ignition. The duration of that potential phase of spin-down, or internal angular-momentum redistribution, is one of the important remaining uncertainties for this scenario. However, a long spin-down phase would be problematic for the ejection of high-velocity helium-burning donor stars. For producing a population of high-velocity He-sdO stars, a natural scenario within the spin-up-spin-down model may be if, in these cases, the rate of ongoing mass accretion becomes insufficient to sustain the rotation needed to avoid explosion (e.g. if the WD reaches a mass at which differential rotation is necessary to delay explosion further).

3. Physical assumptions and methodology

3.1. Numerical method

This study employs the lightweight stellar kinematics code SHyRT (Neunteufel 2020; Neunteufel et al. 2021) in order to model the movement of ejected SN survivors through the Galactic potential, which was assumed to be static and correspond to model 1 put forward by Irrgang et al. (2013) as a revision of Allen & Santillan (1991) in the form

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{b}}(R)=-\frac{M_{\mathrm{b}}}{\sqrt{R^2+b_{\mathrm{b}}^2}}\end{array}$$(1)

for the bulge, where R is the distance from the Galactic centre (note that the gravitational constant G = 1 here)

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{d}}(r,z)=-\frac{M_{\mathrm{d}}}{\sqrt{r^2+(a_{\mathrm{d}}+\sqrt{z^2+b_{\mathrm{d}}^2})}}\end{array}$$(2)

for the disc, where r is the distance from the Galactic centre in the x–y-plane and z is the distance from the x–y-plane and

$$\begin{array}{c}\hfill {\mathrm{\Phi }}_{\mathrm{h}}(R)=\{\begin{array}{cc}-\frac{M_{\mathrm{h}}}{a_{\mathrm{h}}}[\frac{1}{\gamma -1}ln\left(\frac{1+{(R/a_{\mathrm{h}})}^{\gamma -1}}{1+{(\mathrm{\Lambda }/a_{\mathrm{h}})}^{\gamma -1}}\right)-\frac{{(\mathrm{\Lambda }/a_{\mathrm{h}})}^{\gamma -1}}{{(1+\mathrm{\Lambda }/a_{\mathrm{h}})}^{\gamma -1}}]\hfill & \text{if}R<\mathrm{\Lambda }\hfill \\ -\frac{M_{\mathrm{h}}}{R}\frac{{(\mathrm{\Lambda }/a_{\mathrm{h}})}^{\gamma -1}}{{(1+\mathrm{\Lambda }/a_{\mathrm{h}})}^{\gamma -1}}\hfill & \text{if}R\ge \mathrm{\Lambda },\hfill \end{array}\end{array}$$(3)

with γ = 2 and the other parameters given in Table 1.

This potential is axisymmetric and neglects Galactic substructure such as spiral arms and the bar in the disc as well as cold dark matter substructure in the halo. This choice may affect the ratio of bound to unbound objects on long timescales, but the effect is deemed unimportant for our simulation time limit of 300 Myr. Specifically, we analyse six different simulation runs, one each for one assumed terminal WD mass. Parameters are given in Table 2. Note that run M14 was not performed specifically for this study, but Neunteufel et al. (2021), but we present a more detailed and extended analysis here. We draw attention to the fact that the upper boundary for the mass of the runaway differs between individual runs. This is a result of the ejection velocity spectra obtained by Neunteufel (2020) only covering the stated ranges. While this introduces some difficulty in comparing the different synthetic populations, we opt to remain within the bounds of the provided spectra, as opposed to presenting an extrapolation or constraining the samples overall. Instead we highlight the differences in sample size and construct comparable sub-samples if they would, quantitatively or qualitatively, affect our results.

Each run contains a total number of N_e individual tracks, with initial positions randomly assigned in the Galactic x–y-plane, but weighted to correspond to the Galactic density profile as defined by model 1 provided by Irrgang et al. (2013) and see Neunteufel et al. (2021) for the initial distribution. The ejection direction is chosen randomly across 4π. The ejection velocity is determined by the randomly chosen mass of the ejected star, randomly falling into the ranges (depending on mass) detailed in Table 3, as provided by Neunteufel (2020).

3.2. Physical assumptions

Each trajectory is followed for a duration of 300 Myr, with one object assumed to be ejected every 300 yr. This ejection rate corresponds to the inferred Galactic SN Ia rate (Cappellaro et al. 1997; Sullivan et al. 2006). Unless explicitly stated otherwise, we characterise each of our synthetic populations as a snapshot of all ejected objects at the end of the simulated 300 Myr time frame. Each object is assigned a predetermined lifetime corresponding to its nuclear timescale at the start of its helium-burning phase. The latter strategy is a result of our attempt to remain agnostic towards the explosion mechanism leading the ejection of the runaway. Different explosion mechanisms lead to a variety of requirements for the terminal accretor mass and the mass accretion rate it is subjected to. It follows that, by accepting the resulting agnostic ejection velocity spectra as input parameters of our study, we lose sensitivity to the remaining lifetime of the ejected star. Keeping this in mind, the nuclear timescale of a He-star of equal mass to the ejected star can serve as a reliable, though overestimated, upper limit. The mass-dependent lifetime assumed in this study is depicted in Fig. 1. The values shown in Fig. 1 were derived from the nuclear timescales of MESA models of the appropriate mass (Paxton et al. 2011, 2013, 2015, 2018, 2019) release 12 778. As such, while all ejection trajectories are followed for the full duration mentioned above, stars exceeding their mass dependent lifetime are flagged appropriately.

Fig. 1. Assumed lifetimes τ of the ejected He-stars. Stars in the shaded region are below the helium-burning limit and are deemed to thermally collapse to form WDs within 10 Myr.

Each terminal donor mass is chosen randomly in a range corresponding to the values given in Table 2. The mass distribution is assumed to be flat. The latter assumption is again motivated by our desire to remain agnostic towards the explosion mechanism, as a non-flat initial mass distribution may be a result of either the explosion mechanism (i.e. certain terminal donor masses are favoured) or of the progenitor-binary’s pre-evolution. Small number statistics in the observational record precludes reconstruction of an empirical initial mass distribution (i.e. the mass distribution of the observed population after ejection), while uncertainties in the theoretical modelling preclude a reliable construction of an initial mass distribution. We do note that population synthesis models predict that the initial mass distribution of the progenitor system (i.e. the mass distribution at the time of formation of the progenitor binary) is not flat (see e.g. Wang et al. 2013). However, as seen in the aforementioned reference, the mass distribution of the donor star strongly depends on assumptions on the common envelope efficiency parameter α, with lower end estimates producing an approximately flat distribution. As the observational record for these objects improves, reconstruction of an empirical initial mass distribution may become feasible, but for the moment we opt to assume a non-biased approach and allow the initial mass distribution to be flat as a baseline. However, we investigate some variation in initial conditions in Sect. 6.

As seen in Table 3, each terminal donor mass is associated with a range of ejection velocities. Each individual ejection velocity is then randomly chosen from this range, again assuming a flat distribution. As shown by Bauer et al. (2019) and, more recently Liu et al. (2021), interaction with the SN ejecta imparts a non-zero kick on the surviving donor star of the order of 200 km s⁻¹. However, as this kick is (in circularised systems, as under consideration here) imparted perpendicular to the star’s orbital motion, the pre-kick orbital motion dominates the post-kick velocity of the runaway (the difference is of the order of 20 km s⁻¹). Since the magnitude of this kick depends on the pre-kick orbital configuration and the particulars of the SN (i.e. orbital separation at the time of the SN and energetics of the explosion), we neglect it in our calculations. The SN progenitor system is deemed to be commoving with the surrounding stellar environment which, in turn, is assumed to orbit the Galactic centre. Therefore, the ejected star is given an additional initial velocity relative to its assumed rotation around the Galaxy as determined by model 1 of Irrgang et al. (2013).

A sample of the resulting trajectories is shown in Fig. 2. It should be noted that, while the number of experiments for each run is the same, the terminal mass of the exploding WD (M_WD, f) is varied. However, since we opted to sample the ejection velocity spectra provided by Neunteufel (2020), which are subject to a M_WD, f-dependent upper limit for the terminal mass of the runaway sdB star, the number of runaway mass samples for each terminal accretor mass is different. These differences between our individual runs, as well as their designations, are summarised in Table 2.

Fig. 2. Illustration of 75 tracks, out of 1×106, chosen at random, from the MWD,f = 1.4 M⊙-population. Coordinates are Galactic Cartesian. Origins, all set at z = 0, are indicated by red circles. Tracks dropping below the Galactic plane are indicated by dashed lines.

In our simulations we neglect of the effects of dynamical friction for runaways crossing the Galactic stellar population. The timescale T_DF for DF to change the speed v of the runaway (mass M) by order itself can be estimated as (Binney & Tremaine 2008, S8.1)

$$\begin{array}{cc}\hfill T_{\mathrm{DF}}& \equiv {\left(\frac{1}{v}\frac{\mathrm{d}v}{\mathrm{d}t}\right)}^{-1}\sim \frac{v^3}{4\pi G^{2}M\rho ln\mathrm{\Lambda }}{[\mathrm{erf}(X)-\frac{2X}{\sqrt{\pi }}exp(-X^2)]}^{-1}\hfill \\ \hfill & \approx \frac{v^3}{4\pi G^{2}M\rho ln\mathrm{\Lambda }}\hfill \\ \hfill & \approx 8.4\times {10}^{17}\mathrm{yr}{\left(\frac{v}{600\mathrm{km}{\mathrm{s}}^{-1}}\right)}^3{\left(\frac{M}{0.4M_{\odot }}\right)}^{-1}{\left(\frac{\rho }{0.1M_{\odot }{\mathrm{pc}}^{-3}}\right)}^{-1}\hfill \\ \hfill & \times {\left(\frac{log\mathrm{\Lambda }}{27}\right)}^{-1},\hfill \end{array}$$(4)

with $X\equiv v/(\sqrt{2}\sigma )$ and σ as the velocity dispersion, ρ is the stellar density, and Λ is the Coulomb factor, estimated as Λ ∼ b_max/b₉₀ ∼ b_maxv²/[G(M + m_⋆)] with b_max being the maximum impact parameter and m_⋆ the average stellar mass. Assuming M = 0.4 M_⊙, m_⋆ = 1  M_⊙, v = 600 km s⁻¹, which is a reasonable lower limit for stars ejected to unbound orbits, and b_max = 10 kpc, logΛ ∼ 27. The second line in Eq. (4) applies since X ∼ 11 assuming σ = 40 km s⁻¹ for the Solar neighbourhood (e.g. Binney & Tremaine 2008), implying that the factor in square brackets is extremely close to unity. For these assumed values and setting ρ = 0.1  M_⊙ pc⁻³ for the Solar neighbourhood (e.g. Binney & Tremaine 2008), Eq. (4) yields that T_DF exceeds the Hubble time by many orders of magnitude, thus showing that dynamical friction for unbound objects is unimportant.

4. Results – kinematic overview

We performed 6×10⁶ experiments as described above. However, in order to keep the length of this article manageable, we opt to focus our analysis on subsets of our data as required by the current topic of discussion. Specifically, we discuss general results (i.e. results easily transferred to the rest of the set) by using the M14 set as a baseline. Where necessary we expand to the M10, M12 and M15 sets as a comparison. The M11 and M13 sets are not discussed in detail and only mentioned where required. For each synthetic population, we assume one ejection (i.e. one SN producing a runaway) every 300 years, in line with estimates of the inferred Galactic SN Ia rate (Cappellaro et al. 1997; Sullivan et al. 2006). Each trajectory is then randomly assigned an initial ejection time and its position at the end of a 300 Myr period is recorded. The final sample for each synthetic population is composed of the position and velocities, as well as the mass and state (pre-WD or WD) for each individual trajectory.

4.1. Trajectories

The velocity of a star ejected in a SD-He SN is determined mainly by its terminal orbital velocity, to which is added the local orbital velocity of the commoving stellar environment. In the following, we discuss the correlation of general properties of population orbits as measured at a given location.

Here, as a practical choice, we measure the velocity field at two locations. Once at the position of Earth with respect to the Galactic disc, and once at the location of US 708 in the same frame. A full discussion of the velocity field is possible given our data set, but is deemed beyond the intended scope of this study (and deemed more worthwhile if the observational record improves in the future). For this investigation we define a radius of 1 kpc around the locale of observation as “local space”. Figure 3 shows all trajectories (independent of ejection time) in the M14 data set which, at some point during the simulated time interval, cross into Earth’s local space (though trajectories originating within this space are neglected) as well as the ejection locations of each one of these tracks, together with whether they are bound or not. During the following discussion it should be noted that the Galaxy rotates, as seen from above (with Earth located above the Galactic plane), in a clockwise direction (Bhattacharjee et al. 2014). Figure 3A shows a clear preference for unbound stars to leave Earth’s local space (and then the Galaxy) in a prograde fashion with respect to Galactic rotation, while bound objects have a preference for retrograde motion.

Fig. 3. Trajectories of all stars crossing into Earth local space from the outside. Unbound trajectories are indicated as red lines, bound trajectories as blue lines. Panels A and B: galactic Cartesian coordinates in the x–y-plane, with panel A showing the full trajectories and panel B the respective point of origin. Panel C: full trajectories in the x − z-plane. The positions of Earth and the Galactic centre are as indicated. Note that, as seen from above the Galactic plane (with Earth located above the plane), the Galaxy rotates in a clockwise fashion.

Figure 3B shows the ejection locations for the same tracks. As can be seen, trajectories crossing into local space are dominated by stars ejected from the vicinity of Earth, with tracks originating from prograde locations dominated by bound stars and tracks originating from retrograde locations dominated by unbound stars. The Galactic centre, where, due to the higher assumed density of progenitor binaries, a higher number of stars are assumed to be ejected, plays a secondary role compared to Earth’s immediate vicinity. There is further a preference for tracks originating from locations at shorter distance from the Galactic centre. This is in line with the expectation that higher density regions dominate. Figure 3C shows the Galaxy in a the x − z-plane. As seen, the tracks leaving the Galaxy in a positive z-direction dominate. This is a consequence of or initial assumption of all ejection locations being uniformly located at z_Gal = 0 kpc but Earth being located at about z_Gal = 0.4 kpc. If ejection locations were (more realistically) scattered throughout the volume of the disc, this preference would likely diminish, but not vanish.

The correlation between kinematics and ejection location is further exemplified by Fig. 4. This figure shows the ejection location of each star in the Galactic plane, given in terms of angle centred on the Galactic centre with Earth located at 180 deg. This means that, the Galaxy seen from above, the left half of the plot shows the lower (retrograde) half of the Galactic disc while the right half of the plot indicates the upper (prograde) half. The upper half of the plot shows ejection direction aligned retrograde to Galactic rotation, the lower half prograde. The y-axis also provides a notion of the distance of the ejection location from the Galactic centre. This figure is thus divided into four quadrants. Unbound trajectories crossing into Earth’s local space are clearly correlated with retrograde ejection locations and prograde ejection directions while bound trajectories tend (though not quite as strongly correlated) to originate from prograde ejection locations while having their ejection directions aligned in a retrograde fashion. The retrograde-retrograde quadrant is unsurprisingly dominated by bound trajectories, which orbit the Galaxy once (in a retrograde direction) before crossing into Earth’s local space¹. As we do not take the movement of Earth into account in this investigation, the observer would have moved away from the sample location by the time stars on the retrograde-retrograde quadrant orbit the Galaxy. The lower right, prograde-prograde, quadrant is largely empty. The conclusion from this is that stars ejected prograde from a prograde location do not tend to cross into Earth’s local space, except for a very small part of the parameter space where stars are ejected largely radially. The absence of any unbound trajectories in this part of the parameter space is not surprising, as these stars will simply be ejected. More surprising is the absence of any bound trajectories, which can be interpreted as Galactic rotation, all else being the same, being a decisive factor in whether an ejected star is bound or unbound. Further, stars on bound trajectories of this type will likely require more than the interval 300 Myr taken into account in our simulation to cross back into Earth local space. We further find that a star approaching Earth from a prograde position on a retrograde orbit is likely bound and vice versa.

Fig. 4. Chirality of the ejected stars crossing into Earth local space in the M14 data set with respect to Galactic rotation. On the x-axis (location), prograde and retrograde refer to the ejection location of the star being located either prograde or retrograde in the Galactic disc with respect to Earth. On the y-axis (direction), prograde and retrograde refer to the initial velocity vector of the ejected star being oriented either prograde or retrograde with respect to Galactic rotation. The position of Earth is as indicated. Blue dots imply bound stars, red unbound stars.

Similarly, Fig. 5 shows trajectories intersecting a 1 kpc volume surrounding the current position of US 708 (i.e. US 708 local space). As seen, there is a clear preference for trajectories intersecting US 708 local space to be unbound rather than bound. This is equivalent to the statement that a star is more likely to reach a position such as that of US 708 if it is unbound rather than bound. As in Earth local space, unbound trajectories tend to leave the Galaxy in a prograde direction. There is a clear preference for stars to reach the position of US 708 if they are ejected perpendicular to the Galactic plane (note that this is not the case of US 708 itself). As further seen, the ejection location of US 708 is not clearly favoured over locations below (on the positive z-axis) its current position. As for Earth local space, there is a clear preference for ejections towards the positive z-axis to reach US 708 local space.

Fig. 5. Same as Fig. 3 for the US 708 local space. Here the inferred current positions of US 708 and its ejection location are indicated as well.

4.2. Number density

The expected spacial distribution of the runaway population is modelled according to the precepts set out in Sect. 3. The population overview shown here represents a snapshot of the population at the end of the 300 Myr simulation timeframe.

In Fig. 6 we construct a two-dimensional histogram along Galactic cylindrical coordinates for the M10, M12, M14 and M15 populations. This illustrates a difference in distribution between the progenitor population which, (if initially distributed in the volume of the disc as opposed to only in the plane) would follow the black isopycnics. For the ejected population, the local space-number density is depicted in colour scale with each change in colour indicating a change in magnitude, and each black line indicating a drop in magnitude of baryonic matter density, starting at the Galactic centre. It is evident that the density gradient of the ejected population is shallower than for the progenitor population. At this point it is evident that the structure of the ejected population depends on the terminal accretor mass of the progenitor binary, with the lowest terminal accretor mass, in the M10 population, correlated with a more spherical space distribution, while higher terminal accretor masses would produce a flat, cylindrical space distribution. In general, all other parameters being the same, lower terminal accretor masses lead to a larger number of runaways remaining in the vicinity of the Galaxy. Distinctively, synthetic populations with higher terminal accretor masses (M12, M14 and M15) exhibit a roughly cone-shaped volume directly above and below the Galactic bulge largely devoid of stars. While also present, though less pronounced, in the M10 population, this feature is a result of Galactic rotation leading to radial ejection directions being favoured over ejection directly perpendicular to the Galactic plane.

Fig. 6. Expected number density distribution per kpc−3 in cylindrical galactocentric coordinates for different terminal accretor masses with MWD,f = 1.0 M⊙, MWD,f = 1.2 M⊙, MWD,f = 1.4 M⊙ and MWD,f = 1.5 M⊙ depicted in panels A–D respectively. Here, the colour scale indicates the calculated number density. The current positions of Earth and US 708 are indicated. Black lines indicate surfaces of constant baryonic matter density, logarithmic in base 10, as proposed by Irrgang et al. (2013), model 1. White areas do not contain ejected stars in this realisation of our model.

Comparison of the number density of objects in the vicinity of Earth and in the vicinity of US 708, as summarised in Table 4 and Fig. 7, reveals, apart from a higher expected number density in the vicinity of Earth, an approximately constant relative density between the two locations, see green dashed line and error bars in Fig. 7². One of the possible reasons that US 708 was discovered before any comparable objects closer to Earth may be (apart from the conclusion that they do not exist) that US 708-like objects could be difficult to disentangle from a denser stellar environment. Another possible reason, namely that US 708 (like all other candidates) is extreme in at least one observable, will be discussed in Sect. 8.

Fig. 7. Comparative number densities in the vicinity of Earth and US 708 as described in Table 4. Error bars indicate one σ.

5. Qualitative properties of the predicted population

5.1. Defining criteria

Taking into account the theoretical stellar properties of population members, observational evidence and kinematic properties as discussed here, a star strictly qualifies as a candidate member of this population by fulfilling the following criteria:

1. One of:

   • (a)

     Strong He-enrichment with log n_H/n_He ≲ −2, ideally spectroscopically similar to US 708.

   • (b)

     Spectra compatible with a proto-WD.

   • (c)

     Currently a WD.

2. Either of:

   • (a)

     Currently located in the disc.

   • (b)

     If part of the UI or BI subpopulations (see Sect. 5.2): Time of flight to the most recent disc intersection compatible with typical sdB-lifetimes.

3. Ejection velocities at the most recent disc crossing in the range 400−1150 km s⁻¹, depending on mass (if determined).

4. Likely mass < 1.0 M_⊙.

5. Time since ejection < 300 Myr.

We note that item 5 is only necessary as a consequence of the constraints adopted for our simulation. Bound WDs, i.e. descendants of bound non-degenerate runaways, may have been ejected earlier than in the given time frame. Specifically, an inferred ejection velocity which includes the Galactic centre does not disqualify candidates. Distinction from HVS produced by interaction with the SMBH (Hills 1988, and similar) at the Galactic centre should chiefly rely on velocity and structure arguments, not the inferred ejection location, unless the trajectory can be traced back there directly (as in Koposov et al. 2020). Further, being unbound from the Galactic potential should also not be used as a criterion, as the question of boundness relies heavily on the intrinsically random parameter of the ejection direction. Bound and unbound objects instead form their own subpopulations.

5.2. Subdivisions

5.2.1. Boundness

As shown earlier in Sect. 4 and later in Fig. 8C, in the benchmark model, runaways of masses M_*, run < 0.3 M_⊙ will not produce bound objects, while all other runaway masses may result in both bound and unbound objects. Since in our model further interactions with the remaining stellar population are neglected, an initially bound object will remain bound for the duration of our reference time. While, as stated, the question of an observed object being bound is, due to uncertainties in astrometry and insufficiently precise knowledge of the Galactic potential, not easily answered, boundness serves as a convenient division of the population. As such we designate objects bound to the Galactic potential (these would be called halo stars or extreme halo stars observationally) as “B-type”. Conversely, unbound objects are designated “U”. As will be shown later, terminal accretor masses M_WD, f < 1.4 M_⊙ may produce B-type objects at all considered runaway masses. It should be mentioned that, if the Galaxy were to considered as a more realistic N-particle system instead of a static potential, as in this study, multi-body interaction post-ejection may lead to B-type objects becoming U-type and vice versa. However, since the size of the total population (as well as important parameters such as its mass spectrum) is currently not discernable observationally, this possibility is deemed unimportant for this study.

Fig. 8. Spectral distribution of subpopulation membership over terminal donor mass for the M10, M12, M14 and M15 populations for the total set (not merely the fiducial volume). Total numbers per mass at measured against the left hand y-axis. The black line indicates the BI/UI and is measured against the right hand y-axis.

5.2.2. Evolution

As shown by Neunteufel (2020) the descendant of the donor star of the progenitor binary is unlikely to have the same appearance as a WD. Instead, if its mass at the time of ejection is larger than 0.3 M_⊙, it is likely to appear observationally as a He-star (likely very similar to a He-sdO, like US 708). While some introduction of additional entropy into the envelope of the ejected object cannot be ruled out, the star is expected to return to its previous state after approximately one Kelvin-Helmholtz timescale (Liu et al. 2015, 2021; Bauer et al. 2019). The star, then a runaway, will continue its evolution as a He-star, eventually becoming a WD. If the mass of the ejected component prior to ejection decreases to below 0.3 M_⊙ the star becomes unable to sustain He burning in its core, and its state at this point is more precisely described as a proto-WD. However thermal contraction during this stage proceeds slow enough that, under the influence of GWR emission, the system never detaches, which leads to the conclusion that, at ejection, the former donor star must still be in the state of a not fully degenerate proto-WD. Further, once the proto-WD reaches a mass ∼0.2 M_⊙, angular momentum transfer to the binary companion dominates of angular momentum loss through GWR, leading to an increase of the orbital period and a decrease of the ejection velocity Tutukov & Yungelson (1979), Yungelson (2008), Neunteufel (2020) For these objects, as shown by Bauer et al. (2019), Liu et al. (2021), interaction with the SN ejecta will lead to significant loss of material as well as an increase in entropy of the remnant, leading to an increase in radius as well as luminosity. The excess energy is then radiated away and the star contracts, analogously to the higher-mass case, thermally. However, unlike in the higher-mass case, the proto-WD does not return to its previous state. Instead it simply collapses to form a proper WD. As such, at ejection, all population members will be either core He-burning stars or proto-WDs, evolving to form degenerate objects independently (pre-WDs). Therefore the population naturally consists of pre-WDs (I-type stars) and WDs (II-type-stars), giving rise to the second subdivision. Population member therefore evolve naturally from being I-type objects to the state of II-type objects.

5.2.3. Expected subpopulations

The two criteria established above give rise to four distinct subpopulations: firstly, bound pre-WDs, designated BI; secondly, bound WDs, designated BII; thirdly, unbound pre-WDs, designated UI; and lastly, unbound WDs, designated UII. As discussed, in our model, objects are assumed to evolve from subpopulation BI to BII and UI to UII, but not the reverse, but not from BI to UI or the reverse or BII to UII or the reverse. However, as stated, since stars of masses < 0.3 M_⊙ are structurally different, and evolve differently, to those above this limit, we distinguish high mass (> 0.3 M_⊙) UI and BI objects and low mass (< 0.3 M_⊙) BI and UI objects.

6. Results – quantitative properties

This section is dedicated to a quantitative overview of the predicted population. We further discuss some observational properties, as well as some bulk properties of our sample. As in Sect. 4.2, our analysis in this section is based on a snapshot of our population taken after 300 Myr.

6.1. Fiducial volume

Naked helium stars of even moderate masses can be quite luminous (Heber 2009, 2016) and are thus detectable out to considerable distances (∼8.4 kpc in the case of US 708). However, their luminosity is similar to intermediate mass main sequence stars, making helium stars in regions of significant reddening, such as the disc (Geier et al. 2019), liable to be misclassified as MS stars (Irrgang et al. 2019). Further, the part of Hertzsprung-Russell diagram populated by sdO-B stars is shared with other objects, such as hot WD stars. Thus, spectral analysis is the only way to reliably classify these objects. Accordingly of the ∼40 000 candidates reported by Geier et al. (2019), only about 10% have been spectrally classified. See also Neunteufel et al. (2021) for a discussion. In order to facilitate comparison of our predictions with the observational record, we reject members of our synthetic population lying outside a pre-determined region deemed inaccessible to observation. Specifically, we assume that the likelihood of (eventually) detecting a member of our predicted population is unity inside our fiducial volume and zero everywhere else. We define this fiducial volume thus to consist of a 2 kpc sphere around the position of Earth and, additionally two cones, their axes of symmetry perpendicular to the Galactic plan and their vertices at Earth, one pointing up and the other down (with Earth above the plane). The angle between the surface of both the upward and downward pointing cone and the plane is α = 10°, both cones extending out to a radial distance of 25 kpc from Earth. The limitation on the 2 kpc sphere is motivated by the approximate limit for obtaining reliable Gaia-parallaxes for very blue objects such as the ones under consideration here. The value for α is approximate for limits of negligible reddening as indicated by the results of Green et al. (2018).

The shape of our fiducial volume is illustrated in Fig. 9. We note that, as a consequence, any computation of the observationally supported production rate based on the arguments presented here should strictly reject any new observation outside of the fiducial volume as defined above. We further note that the shape of the fiducial volume is only an approximation. The geometry of Galactic extinction and reddening at high Galactic latitudes is non-trivial (Green et al. 2018; Sun et al. 2021).

Fig. 9. Illustration of the fiducial volume chosen for the calculation of the expected detection and event rates. Only objects located in the red shaded area are taken into account. Specifically, all objects in a distance of D < 2 kpc are assumed to be visible. Any objects at a distance D > 25 kpc are assumed to be undetectable. Objects located in the disc are discounted by the assumption that objects with a Galactic latitude α < b < −α are invisible, with α = 10°. Earth’s position is as indicated. Illustration not to scale.

6.2. Cumulative distributions

Figure 10C shows the sum of cumulative distributions of the UI, BI and BII subpopulations (the UII subpopulation is likely undetectable by virtue of being composed of fast, low mass WDs at large distance) in the M14 population with respect to distance from Earth. As can be seen, the BII subpopulation constitutes of the order of 10% of the total population. The BI subpopulation constitutes of the order of 40% at distances smaller than 15 kpc, decreasing to about 30% at 25 kpc. Within a volume extending to the distance of US 708 (itself a member of the UI subpopulation) about half of all objects are part of the UI subpopulation. The distribution resulting from the M10 population is shown in Fig. 10A. As can be seen, about 77% (as opposed to about 40% in the M14 sample) of the population is constituted by the B-subsample, with about 40% being composed of the BII subpopulation and the remaining 37% of the BI subpopulation. We note that this, to some extent, is a result of pollution by runaways in the 0.8−1.0 M_⊙-range, which are not present in the M14 sample. These objects account for a pollution of about 10% in the BI subpopulation and about 28% in the BII subpopulation, illustrating the importance of the terminal donor mass for the structure of the resulting population. discounting the higher mass pollution, the BI subpopulation is about 50% of the total population and the BII subpopulation about 20%, the UI subpopulation accounting for the rest. It is clear, however, that the ratio of members of the different subpopulations, UI to BI+BII, is strongly influenced by the terminal accretor mass. The distribution for the M12 population, as shown in Fig. 10B, falls between the M10 and the M14 set (as would naively be expected). Here, the UI-subpopulation constitutes about 66% of the total population (including the 0.8−1.0 M_⊙-admixture). This provides further evidence of the terminal accretor mass being the decisive factor in the determination of boundness to unboundness. For the M15 population, shown in Fig. 10D, a clear preference for the UI population when compared to the other synthetic populations (only comparing the mass ranges present in all populations) becomes evident. As will be discussed in Sect. 6.3, the B–U ratio may constitute a viable criterion to distinguish terminal accretor masses. A clear conclusion of the discussion in this section should be the realisation that the relative number of members of the UI-subpopulation compared to the B-subpopulation is strongly influenced by the terminal accretor mass of the generating SN. Further, if a significant number of runaways have a mass exceeding 0.7 M_⊙, then a population of high-velocity WDs in the Galactic halo (i.e. the BII-subpopulation) is a necessary consequence. We mention in passing that this prediction is similar to that by Hansen (2003), who considered a slightly different channel, requiring a donor in the mass range 1.3−3.3 M_⊙, leading to lower ejection velocities incapable of producing HVS.

Fig. 10. Cumulative distance distribution of runaway stars normalised to the total number of runaways in a non-degenerate state, Ntot, n − d for a ejection frequency of fej  =  1/300 yr−1. The blue and orange curves add up to unity, with blue indicating the relative abundance of the UI subpopulation and orange represents the BI subpopulation. Red indicates the BII subpopulation normalised to the UI+UII subpopulation. Dark orange and dark red indicate the admixture of objects in the mass range 0.7−0.8 M⊙. Panels A–D: M10, M12, M14 and M15 populations respectively.

6.3. The B–U-ratio as a measure of terminal accretor mass

In Fig. 11 we show the relative occurrence rates of individual subpopulations as defined in Sect. 5.2 within the fiducial volume, comparing different synthetic populations, but excluding terminal donor masses of M_⋆, run > 0.7 M_⊙ in order to perform the comparison within a similar parameter space. Since the ratio of the BI to the BII subpopulation is affected by the assumed lifetime of each ejected star, we perform our analysis once under the assumption of our baseline lifetimes and once under the assumption of the baseline being reduced by half. As can be seen, the B/U-ratio is strongly affected by the choice of terminal accretor mass, reducing from about 1.25 at M_⋆, run = 1.0 M_⊙ to about 0.4 at M_⋆, run = 1.5 M_⊙. This indicates that, given a sufficiently large sample, a dominant accretor mass (assuming one dominates) should be discernable in the observational record with the B/U-ratio being a good indicator. Other than the BII/BI-ratio, this ratio is unaffected by any assumption on the runaway lifetime. In the same vein, while the BII/BI-ratio is a function (though a weaker one than the B/U-ratio) of the terminal accretor mass, this ratio is deemed to be a less viable indicator of the terminal accretor mass (though it may be a good indicator of the delay time, as it encodes the remaining lifetime of the donor star, but the necessary analysis is beyond the scope of the paper).

Fig. 11. Relative occurrence of individual subpopulations in all synthetic populations in the fiducial volume depending on different terminal accretor masses. BII/BI(τ/2) was calculated under the assumption of a 50% reduced nuclear burning lifetime of each ejected star.

6.4. Dependence on the ejection rate

In the previous discussion, an ejection rate of f_ej = 1/300 yr⁻¹ was assumed as a baseline, roughly corresponding to the inferred Galactic SN Ia rate (Cappellaro et al. 1997; Sullivan et al. 2006). In this section, we investigate whether changes in this assumption affect the qualitative conclusions drawn from our synthetic population. This investigation is effected by randomly choosing, from our initial synthetic population, a subset of tracks such that the resulting ejection rate is decreased either by a factor of ten, or a factor of one hundred from the baseline assumption. We perform this investigation explicitly for the M14 population. Figure 12 shows the resulting cumulative distribution within the fiducial volume for an assumed ejection rate of f_ej = 1/3 × 10⁻⁴ yr⁻¹. This means that over the same 300 Myr time frame, instead of 10⁶ objects being ejected, only 10⁴ objects are ejected. We note that, despite some statistical fluctuations due to the smaller set of objects, the relative numbers between the different subpopulations are only marginally affected.

Fig. 12. Same as Fig. 10 with fej = 1/3 × 10−4 yr−1.

We performed a comparison with a similar calculation assuming f_ej = 1/3 × 10⁻³ yr⁻¹. The number of objects contained in the fiducial volume, N(D = 25 kpc), decreases by a factor of ∼10.02 between f_ej = 1/3 × 10⁻² yr⁻¹ and f_ej = 1/3 × 10⁻³ yr⁻¹, by a factor of ∼9.61 between f_ej = 1/3 × 10⁻³ yr⁻¹ and f_ej = 1/3 × 10⁻⁴ yr⁻¹ and, accordingly, by a factor of ∼96.27 between f_ej = 1/3 × 10⁻² yr⁻¹ and f_ej = 1/3 × 10⁻⁴ yr⁻¹. It is therefore reasonable to assume that the number of objects in the fiducial volume decreases linearly depending on the ejection rate, while the relative numbers in each subpopulation are, to first approximation, constant. We assume such a linear dependence, as well as constancy of the relative numbers, in the following discussion. At this point we note that SN Ia rates have been argued to be inversely correlated with the delay time via $\sim \frac{1}{t_{\mathrm{delay}}}$ (Ruiter et al. 2009; Maoz & Mannucci 2012), which implies a time dependence of the ejection rate which is at odds with our assumption of a constant ejection rate. However, we note that, as t_delay = t₀ + t_current, this effect would lead to a decrease of the ejection rate of the order of 2−3% over the 300 Myr interval covered by our simulations. A decrease of this nature has a negligible impact on our results in every meaningful interpretation.

6.5. Dependence on the terminal donor mass

In this section we investigate the contributions of different assumed terminal donor masses to the total and subpopulations. Figure 8 shows the contributions of each mass bin to the different subpopulations as defined in Sect. 5.2 in the M10, M12, M14 and M15 populations (subplots). We note that the flat initial mass distributions. Here, the entire population, independent of location in the vicinity of the Galaxy and ejection timing, is taken into account. The shape of the UI+BI subpopulations relative to the UII+BII subpopulations clearly reflects the assumed lifetimes (see Fig. 1).

Consequently, as expected, the relative contributions of the BII subpopulation increases with the terminal donor mass, provided the terminal donor mass is greater than 0.35 M_⊙. The expectation that there is no contribution of the BII subpopulation in the mass range 0.35−0.4 M_⊙ is a consequence of the assumed lifetime of these objects exceeding the simulation time of 300 Myr. However, this assumption represents an upper limit. In reality, a contribution of the BII subpopulation should be expected in this mass range (as supported by our τ/2 test case). A striking result here is that we expect no contribution of the bound subpopulation (both BI and BII) in the 0.2−0.3 M_⊙ mass range at any terminal accretor mass and only a contribution at terminal donor masses of 0.35 M_⊙. This indicates that the Galactic potential should not be sufficient to bind runaways at this mass range if ejected via the SD-He channel. The only observed candidate object likely falling into this mass range (LP 398-9, see Sect. 7), is very likely unbound as well. If bound objects at the lower end of the mass range (< 0.35 M_⊙) are found, this would be an indication either that the production mechanism differs from the one under consideration here or that the Galactic potential has been underestimated. This conclusion is necessarily independent of the assumed lifetime as defined in Fig. 1, as these objects are proto-WDs and are thus observationally distinct from putative higher mass objects, first presenting as proto-WDs, then as ELM-WDs. Figure 8 further shows the dependence of the ratio of the BI subpopulation to the UI subpopulation. This ratio decreases (at M_⋆, run = 0.7 M_⊙) from ∼0.9 in the M10 population to ∼0.2 in the M15 population, thus showing a clear dependence on the terminal accretor mass. This ratio is also largely independent on the assumed runaway lifetimes, as these only affect the ratio of the UI+BI to UII+BII subpopulations. The dependence of the ratio of the UI+BI to UII+BII subpopulations on the lifetimes can be estimated the following way. The total number of objects in the population is simply the ejection rate times the simulation interval t_sim

$$\begin{array}{c}\hfill N=f_{\mathrm{ej}}·t_{\mathrm{sim}}.\end{array}$$(5)

Then the number of objects in a certain terminal donor mass bin is

$$\begin{array}{c}\hfill N(M_{\star ,\mathrm{run}})=f_{\mathrm{ej}}(M_{\star ,\mathrm{run}})·t_{\mathrm{sim}},\end{array}$$(6)

where, due to the assumed flat mass distribution f_ej(M_⋆, run) = f_ej/n_bins. Then the number of objects in the UI+BI subpopulation is, with the assumed lifetime τ(M_⋆, run), considering a constant f_ej

$$\begin{array}{c}\hfill N(M_{\star ,\mathrm{run}},I)=f_{\mathrm{ej}}(M_{\star ,\mathrm{run}})·\tau (M_{\star ,\mathrm{run}})\end{array}$$(7)

and the number of the UII+BII subpopulation given by N(M_⋆, run, II) = f_ej(M_⋆, run)⋅(t_sim − τ(M_⋆, run)), the ratio is

$$\begin{array}{c}\hfill \frac{N(M_{\star ,\mathrm{run}},I)}{N(M_{\star ,\mathrm{run}},II)}=\frac{\tau (M_{\star ,\mathrm{run}})}{t_{\mathrm{sim}}-\tau (M_{\star ,\mathrm{run}})}.\end{array}$$(8)

As such, while inaccuracies in the assumed lifetimes are suppressed, the BI/UI ratio of the total population is a yardstick for the determination of the terminal accretor mass. As will be shown later, the BI/UI ratio in the fiducial volume, on the other hand, does depend on the assumed lifetime.

As the entirety of the population is spacially inaccessible to observation (as discussed in Sect. 6.1), the population distribution needs to be studied only within the fiducial volume. We show the resulting population in Fig. 13. As seen, the clear dependence of the BI/UI ratio seen in the full spacial sample is replicated in the fiducial volume subsample. A number of features should be noted here: Firstly, the population size increases with increasing M_⋆, run. This is a consequence of the slower ejection velocity of higher mass donor stars. Thus, higher mass donors tend to remain within the fiducial volume for a longer time. Secondly, the BI/UI ratio strongly increases with M_⋆, run. This is, again, a result of the slower ejection velocity of higher mass donors. However, the ratios increase more strongly with decreasing M_WD, f. For instance, not surprisingly, the slower the ejection velocity, the more likely an ejected star is bound. Thirdly, there is a maximum in the BI+UI distribution (not mirrored in the total population) which depends on the terminal accretor mass. While the increase of the population size is driven by the same effect as the continuous increase of the total population (as discussed in item 1), the decrease at higher mass (leading to the maximum) is a result of the decreasing lifetime of the ejected stars as a result of their higher mass. The resulting maximum in the BI+UI subpopulation depends on the terminal donor mass and leads to a maximum at M_⋆, run = 0.55 M_⊙ at M_WD, f = 1.0 M_⊙, M_⋆, run = 0.6 M_⊙ at M_WD, f = 1.2 M_⊙, M_⋆, run = 0.7 M_⊙ at M_WD, f = 1.4 M_⊙; and M_⋆, run = 0.7 M_⊙ at M_WD, f = 1.5 M_⊙.

Fig. 13. Like Fig. 8 but restricted to the fiducial volume. Panel A: the BI/UI ratio was cropped in order to preserve legibility. The missing value is NBI/NUI = 23 at M⋆, run = 1.00 M⊙.

For clarity, Fig. 14 shows only the BI/UI ratios as presented in Fig. 13. As can be seen, the BI/UI-ratio is largely dependent on the terminal accretor mass (as designated by the synthetic population in this figure). The increase of the ratio with decreasing M_⋆, run should be noted.

Fig. 14. Direct comparison of the BI/UI-ratio as presented in Fig. 13.

In Fig. 15 we investigate the effect of a change of the assumed lifetime of the ejected star (τ) by setting τ → τ/2. As this affects the higher mass stars more than the lower mass ones, this change has the consequence of suppressing the UI+BI population vs. the UII+BII population at higher masses, leading to a shift of the peak in the distribution compared to the baseline assumption towards lower masses. However, since all of our trajectories originate in the disc, a star necessarily has to be able to reach fiducial volume within its remaining lifetime in order for it to bee seen as a UI+BI subpopulation member. As such, the BI/UI ratio is affected as well. While, naively, one would expect that this would lead to a decrease of the BI/UI-ratio, we instead, at lower terminal accretor masses (M_⋆, run < 1.3 M_⊙), find an increase. We attribute this to stars ejected in the vicinity of Earth having a higher likelihood of reaching the fiducial volume and stars ejected at larger distances from Earth not being sufficiently fast, even if unbound, to reach the fiducial volume.

Fig. 15. Like Fig. 13 but with assumed lifetime τ = τ/2. Panel A: the BI/UI ratio was cropped in order to preserve legibility. The missing value is NBI/NUI = 55 at M⋆, run = 1.00 M⊙.

We compare the BI/UI-ratios for Fig. 15 with those presented in Fig. 13 in Fig. 16. As seen, the effect on the BI/UI-ratios is largest at lower M_WD, f, becoming negligible at M_WD, f = 1.4 M_⊙ (i.e. in the M14 population).

Fig. 16. Direct comparison of the BI/UI-ratio as presented in Fig. 15. Grey lines indicate the data presented in Fig. 14.

While we assumed a flat initial distribution for M_⋆, run for generality, we investigate the effect of a sloped distribution in Fig. 17 for the M14 population. Here, instead of a flat distribution, we arbitrarily impose a linearly decreasing weighting function (Eq. (9)), which decreases the relative abundance of objects relative to M_⋆, run = 0.2 M_⊙ to zero at M_⋆, run = 1.05 M_⊙ (we note, however, that the most massive runaway star in this synthetic population has a mass of 0.8 M_⊙).

$$\begin{array}{c}\hfill N_{\mathrm{ej}}(M_{\star ,\mathrm{run}})=N_{\mathrm{ej}}(M_{\star ,\mathrm{run}}=0.2M_{\odot })\ast (\frac{-M_{\star ,\mathrm{run}}}{0.85M_{\odot }}+1+\frac{0.2}{0.85}).\end{array}$$(9)

Fig. 17. Like Fig. 13C, but assuming a weighting function according to Eq. (9).

As would be expected, the weighting has the effect of suppressing the population at higher values of M_⋆, run, leading to a shift of the BI+UI subpopulation maximum to lower values. However, we note that the BI/UI ratio is unaffected by this change. The conclusion here is that, should the initial mass function differ from our assumptions (which is likely), and given that the population results from a single terminal accretor mass (which is currently unknown) the BI/UI-ratio would be independent of the relative terminal mass distribution of the donor stars.

6.6. SNR association

With each ejection of a HVS in the SD-He channel, a supernova remnant (SNR) is necessarily created. As such, the observation of a comparable number of SNRs may provide some insight into the viability of our predictions. In Table 5 we compare the number of runaways in each synthetic population with the number of associated SNRs. This is done by counting all runaways within the fiducial volume with a time of flight since ejection short enough for an associated SNR to not have dissipated. The lifetime of a SNR is subject to an number of uncertainties, such as the density of the local interstellar medium, ejecta mass and a number of other parameters. Estimates in literature (e.g. Frail et al. 1994) are of the order of 60 kyr. A back of the envelope calculation reveals that, therefore, over a 300 Myr timespan, about 0.02% of all stars should still have an associated non-dissipated SNR. As such, we assume an upper limit for the average SNR lifetime of 200 kyr. We point out, however, that the expansion velocity of SNRs is generally faster than the highest ejection velocity attainable in this ejection scenario. Therefore, any HVS associated with an SNR should still be physically located inside its SNR except at late stages of SNR evolution. Surveys aimed at finding the ejected donor star in putative SN Ia remnants have generally failed to produce results (Kerzendorf et al. 2012, 2014, 2018a,b), though Ruiz-Lapuente et al. (2004) claimed to have found a candidate for Tycho’s SNR. A notable exception to this is one of the candidate objects for the population under consideration here, namely LP 398-9 (see more detailed discussion in Sect. 7), which has been claimed to be associated with the SNR G70.0–21.5 (Shen et al. 2018), discovered by Fesen et al. (2015).

As can be seen in Table 5, the number of objects in the fiducial volume dominates over the number of associated SNRs by at least a factor of 10⁴, which is on the same order of magnitude as the estimate presented earlier. Remembering that the assumed ejection frequency of f_ej = 1/3 × 10⁻² yr⁻¹ is likely too high, we conclude that non-observation of associated SNRs is the expected outcome even for ideal circumstances. If SNR G70.0–21.5 is actually associated with LP 398-9, this observation should be viewed as coincidental.

7. Candidate members of the population

The observational record of objects fitting into this channel is, at the time of writing, extraordinarily small, with only a single object having been determined to fit with any reliable degree of certainty. However, using the predictions presented here and by Neunteufel et al. (2021), it is possible to determine at least two additional candidate objects for this population. In this section, we briefly discuss the extant observational record. Keeping in mind the arguments presented in Sect. 3, we strive to be inclusive, rather than exclusive, in the construction of this sample. As such, we deem it useful to discuss the entirety of the current sample of candidate objects in detail. In the following, we follow up some astrometric properties of the objects under discussion using the Gaia DR2 and EDR3 data sets (Gaia Collaboration 2016, 2018, 2021).

7.1. US 708

First identified as a faint blue object showing an ultraviolet excess at high Galactic latitude by Usher et al. (1982), then rediscovered as SDSS J093320.86+441705.4 and described by Hirsch et al. (2005), US 708 was further determined it to possess an exceptionally high radial velocity. Justham et al. (2009) argued that a natural explanation for the properties of US 708 is as a surviving helium donor star from a short-period single-degenerate SN Ia progenitor system and that this seemed more satisfying for US 708 than the previously-suggested scenario of dynamical ejection from the Galactic centre (see also Wang & Han 2009 and Geier et al. 2013). Geier et al. (2015a) brought US 708 to prominence as a compelling candidate survivor from a single-degenerate thermonuclear supernova, largely by demonstrating that the ejection location is traceable to the Galactic disc, with no intersection of the Galactic bulge. This largely excluded ejection from the Galactic centre via interaction with the SMBH located there. Brown et al. (2015) independently published similar conclusions for US 708, although noted that the inferred space velocity for US 708 was so high as to be in tension with published SN ejection models of the time, if US 708 was ejected from the Galactic disc, and so speculated that the pre-SN system which ejected US 708 might itself have been a member of the Galactic halo. Neunteufel (2020) re-examined the past motion of US 708 using Gaia DR2 data, which supported the conclusion that this star was not ejected from the Galactic centre. Neunteufel (2020) further found that the most likely inferred ejection velocity for US 708 was lower by about 100 km s⁻¹ than in Geier et al. (2015a). Nonetheless, as shown in Fig. 18, US 708, if it indeed originates from the channel under consideration here, is an extreme member of its own population when its radial velocity is concerned (however, as stated, the extreme nature of US 708 may well be the reason for its discovery).

Fig. 18. Spectral distribution of predicted radial velocities of the predicted UI+BI subpopulation for different terminal accretor masses within the fiducial volume. Lines of different colours indicate the sum over all runaways larger than the mass indicated in the key. Terminal accretor masses for each panel are Mf, WD = 1.0 M⊙ for panel A, Mf, WD = 1.2 M⊙ for panel B, Mf, WD = 1.4 M⊙ for panel C and Mf, WD = 1.5 M⊙ for panel D.

As shown by Fig. 19, the radial velocity exhibited by US 708 is more compatible with low mass (< 0.3 M_⊙) members of the UII subpopulation. However, taking into account US 708’s position in the HRD (generally associated with He-burning objects such as canonical sdB and sdO stars), as well as the time of flight to the most recent disc crossing (its assumed location of ejection) of ∼12−16 Myr suggest that it is currently He-burning, which suggests both a mass > 0.3 M_⊙ as well as a pre-WD structure. As such, US 708 is the prototype of this population as well as a member of the high-mass UI subpopulation.

Fig. 19. Like Fig. 18, but for the UII+BII subpopulation. Each line indicates the sum over all masses higher than indicated in the key; The secondary peak at about 750 km s−1 indicates the peak of the low mass UII population. The lifetimes of stars in the range 0.3−0.45 M⊙ are too long to contribute to this plot.

If US 708 does have a mass > 0.3 M_⊙, we emphasise that the inferred ejection velocity is only marginally consistent with the systems from Neunteufel (2020) in which the WD explodes below the Chandrasekhar mass. If the true mass of US 708 can be determined, this could therefore provide evidence in favour of super-Chandrasekhar single-degenerate thermonuclear SNe. This would also be indirect evidence in favour of the “spin-up/spin-down” model, although the timescale for spin-down or angular-momentum redistribution in this case should be short compared to the helium-burning lifetime of US 708.

7.2. J2050

SDSS J205030.39-061957.8 (J2050) was first identified as part of the Hyper-MUCHFUSS catalogue (Geier et al. 2015b) and then characterised as as high-velocity bound spectroscopic twin of US 708 by Ziegerer et al. (2017), making it a candidate member of the BI subpopulation. With a time-of-flight of > 113 Myr to a disc-transit in the area x_Gal = 10 kpc < x_{Gal, transit} < x_Gal = 20 kpc and y_Gal = −2 kpc < y_{Gal, transit} < y_Gal = 12 kpc, this object is consistent with a He-sdO of mass < 0.4 M_⊙, which is an expected configuration for members of our population. As concluded by Ziegerer et al. (2017), J2050 orbits the Galaxy in a retrograde direction, which makes it consistent with our prediction for members of the high mass BI subpopulation in Earth’s local volume. However, while the aforementioned criteria make J2050 a promising candidate, one major caveat remains in that the inferred ejection velocity (∼385 ± 82 km s⁻¹, as calculated by Ziegerer et al. 2017) is low for products of the single degenerate He-donor scenario. This may be explained by the ejection also producing a bound remnant (as would be expected in certain production scenarios of Type Iax SNe) or some other unresolved mechanism. However, keeping this caveat in mind, J2050 remains one of the most promising candidates known so far.

7.3. LP 398-9

LP 398-9 was identified in the Luyten-Palomar catalogue and first appeared as a high proper motion star in the NLTT catalogue (Luyten 1995, as NLTT 51732) and again in the LPSM-North catalogue (Lépine & Shara 2005, as LSPM J2138+2522). Most notably, this object was identified by Shen et al. (2018) in the Gaia catalogue, as Gaia DR2 1798008584396457088, Gaia Collaboration (2018) as a candidate survivor of a double degenerate thermonuclear SN³. We note that the authors of that study base their conclusion not least on the argument that an inferred ejection velocity of ∼1200 km s⁻¹, by virtue of being too high, cannot be explained by the single degenerate channel. The contradiction arising by comparison with Neunteufel (2020) is explained by Shen et al. (2018) only taking into account terminal He-star of masses down to 0.5 M_⊙ in the calculation of possible ejection velocities. We note that an ejection velocity of ∼1200 km s⁻¹ is too low to be easily explained by the double degenerate scenario, though terminal donor masses < 0.4 M_⊙ and high effective temperatures ∼100 kK could marginally explain the observed value (as shown very recently by Bauer et al. 2021). This problem is somewhat exacerbated by the ejection velocity derived from the Gaia EDR3 data set being on the order 100 km s⁻¹ lower than the velocity derived from the Gaia DR2 date set (see Bauer et al. 2021). However, the inferred ejection velocity is explained with the assumption of the donor star being a He-star or proto-WD with a terminal mass in the range of 0.2−0.3 M_⊙ and terminal accretor masses approaching (or exceeding) the Chandrasekhar mass, which would make it a candidate member of the UI subpopulation. Further, with the appearance of Gaia EDR3 (Gaia Collaboration 2021), the parallax of this object increased to 1.19 ± 0.07 mas from 1.05 ± 0.11 mas, leading to a corresponding correction of the galactocentric space velocity and a decrease of the inferred ejection velocity. As such, LP 398-9 is a promising candidate member of the low mass UI subpopulation.

We would like to draw attention at this point to the proper motion of LP 398-9 in relation to that of the predicted subpopulation as shown in Fig. 20. Comparison with the expected number of objects at the distance of LP 398-9 leads to the conclusion that, a priori, detection of an object at this distance with high proper motion such as this would be deemed very unlikely. Further, the reliability of the Gaia DR2 data, as utilised by Shen et al. (2018), has been been called into question (Scholz 2018).

Fig. 20. Spectral diagram of the total proper motion of the predicted population for the M14 population. Lines of different colours indicate different distance-limited sub-samples. Dashed lines indicate where the distance-limited sub-sample of the same colour is the same as the next smaller limited sub-sample. Grey areas indicate the proper motion of the three population candidates as labelled. The PMs of US 708 and J2050 relatively similar, so their corresponding lines were unified for legibility. Proper motion errors are of the order of ∼ ± 0.2 mas yr−1 for US 708 and J2050 and ∼ ± 0.1 mas yr−1 for LP 398-9 (Gaia Collaboration 2021).

8. Astrometric properties

In this section we discuss the astrometric properties of the population, comparing them with the properties of the observed candidate objects. We again, like in sections Sects. 4.2 and 6, focus on a snapshot of our populations taken at the end of the 300 Myr simulation time frame.

8.1. Proper motion

In the era of Gaia, the most easily obtainable astrometric parameter is the proper motion. Thus, we investigate the proper motion distribution of our population in Fig. 20 for the M14 population. As would be expected, the proper motion depends strongly on the distance of the object, even if high space velocities are considered. If the entire population in the fiducial volume is considered, the proper motion is expected to be scatter around zero, which is borne out by our sample. However, the maximum of the population moves towards higher values of the proper motion if the sample is restricted to smaller distances. Consequently, the maximum of the sample is located at 0 mas yr⁻¹ in a volume extending to 25 kpc. Restricting the sample to 15 kpc leads to a increase of the maximum of the distribution to ∼4 mas yr⁻¹, restriction to 10 kpc to ∼10 mas yr⁻¹ and restriction to 5 kpc to ∼20 mas yr⁻¹. The distribution approaches zero at ∼100 mas yr⁻¹. We note that the observed objects US 708 and J2050, assuming they are part of this population, both exhibit lower proper motions than would be expected for objects of comparable distances (i.e. the sub-sample restricted to 10 kpc). Most notably, the observed object LP 398-9 has been measured to exhibit a proper motion at the upper end of the distribution. In the case of LP 398-9, this high proper motion is a result of a number of factors. Objects at shorter distances (such as ∼1 kpc for LP 398-9) from the observer will generally exhibit a higher proper motion. Further, high space velocity plays a role, with LP 398-9 being the fastest of these three objects. Together with LP 398-9’s peculiar trajectory with a radial velocity not too different from zero, these factors combine to give LP 398-9 an exceptionally high proper motion. Observation of such an object should therefore, a priori, be judged as extremely unlikely. Detection of members of this population should therefore not rely on a lower proper motion cut-off.

8.2. Radial velocity

The radial velocity is not as easily obtainable as the proper motion, but together with the proper motion, the space velocity, which is a primary observable for this population, is fully defined. In Fig. 18 we show the radial velocity distribution for the entirety of the BI+UI subpopulation within the fiducial volume, depending on M_WD, f and decomposed into different values for M_⋆, run. As can be seen, the maximum of the distribution is affected by the value of M_WD, f, with the maximum increasing from ∼450 km s⁻¹ at M_WD, f = 1.0 M_⊙ to ∼520 km s⁻¹ at M_WD, f = 1.5 M_⊙. The distribution maxima for increasing M_⋆, run is, as would be expected, shifted to lower velocities.

We find that objects with the highest expected ejection velocities (M_⋆, run < 0.3 M_⊙) do, counter-intuitively, not significantly affect to the total distribution in RV, as their short lifetime means that they are removed from the sample before significantly contributing. The importance of the lifetime in this discussion becomes evident when investigating the RV-distribution of the UII+BII subpopulation, as shown in Fig. 19. Here, as a larger fraction of ejected stars remain bound as M_WD, f decreases, the majority of the UII+BII subpopulation is clustered around zero, though exhibiting a peak at positive velocities. This peak is generally located at lower velocities than in the UI+BI subpopulation. As such, we expect this channel to result in a significant population of fast moving WDs in this velocity range. We highlight, however, the marked peak of low mass (M_⋆, run < 0.3 M_⊙) objects at around 1000 km s⁻¹. This peak of fast moving, unbound and degenerate WDs is the final destination in phase space of the lower mass UI+BI objects in Fig. 18.

Comparing with the properties of the observed objects with the predicted populations, we find that US 708, due to its high recession velocity, very likely does not originate from a system with terminal accretor mass M_WD, f < 1.2 M_⊙. However, as seen, both US 708 and J2050 are notable as their distinguishing feature is their radial velocity. Both objects are located either in the negative or positive radial velocity tails of the population (note that in the case of US 708, this was already recognised by Neunteufel et al. 2021). We note that J2050’s current radial velocity is an effect of its highly eccentric orbit with respect to the Galactic centre, not of an exceptionally high ejection velocity. As such, both of these objects are unusual in the sense of their radial motion being aligned with Earth in such a way as to provide either the maximum or minimum RV measurements.

LP 398-9 is likely a member of the low mass (M_⋆, run < 0.3 M_⊙) UI population. As members of the low mass UI population have a very short lifetime, decaying into the UII population quickly, their velocity distribution can be assumed to be identical the UII of the same mass. Thus, LP 398-9 is extraordinary in its radial velocity (∼20 km s⁻¹) when compared to the rest of this population, as can be concluded from Fig. 19. As seen, the low mass UII population (secondary peak in the black line in each panel) peaks at no less than 750 km s⁻¹. We emphasise again that this radial velocity measurement was called into question in literature (Scholz 2018).

8.3. Other population properties

We show the expected distance-space velocity relation of the M14 population in Fig. 21. As can be seen, lower mass objects (< 0.3 M_⊙), while higher velocity, are expected to occur less frequently at higher distances, owing to their short assumed lifetime. In panel A, LP 398-9 can be distinguished as an extremely high space velocity object, but not unreasonable for a member of the low mass UI subpopulation. In general, low mass objects have a clear preference for being unbound over being bound. US 708, with an assumed current mass of ∼0.3 M_⊙ fits very well into the predicted parameters of its own mass range, straddling the line between low mass (< 0.3 M_⊙) and higher mass (> 0.3 M_⊙) objects. Further, the cut-off between bound and unbound objects is clearly distinguishable between panels B and C, with J2050 clearly in the bound region occupied by higher mass objects. We emphasise again that this is not consistent with the expected lifetime of J2050 and a most recent disc crossing time of ∼113 Myr. While the previous three panels contain the entire sample, panel D is the same as panel A but under the assumption of halved lifetimes. However, except for an, expected, reduction of the overall numbers and low and high mass objects being visible for a shorter time, the picture and drawn conclusions remain the same. Repeating this discussion for the other synthetic populations yields no significant changes other than for the M10 population, where US 708 would instead be more likely a member of the low mass UI population (which is in agreement with the conclusions drawn by Neunteufel 2020) and LP 398-9 would be an even more extreme member of the low mass UI population.

Fig. 21. Galactocentric space velocity of I-type subpopulation in the M14 (MWD, f = 1.4 M⊙) population in the fiducial volume over distance D from Earth. Colours indicate the mass of the ejected object. Panel A: entire population. Panels B and C: only contain the UI and BI subpopulation respectively. Panel D: the same as panel A, except assuming τlifetime = 0.5 τlifetime, bm. Grey indicates (with no differentiation according to mass) the II-type subpopulation in panels A and D. Positions of the observed candidates in the parameter space are as indicated. The dashed line indicates D = 25 kpc (i.e. the parameter space left of the dashed line is equivalent to the fiducial volume as defined in Fig. 9).

In order to further study the observational properties of the population we calculate the average, as well as maximum and minimum luminosity and effective temperature of stars in our mass range. This approach is necessary because the exact state of the ejected star at the point of ejection, as well as its further evolutionary lifetime, is obscured to us through our methodology. For this, we use a set of models (based on Neunteufel 2020), calculated with MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) release 12 778. The values given for the temperature, as a function of mass (as this parameter is, neglecting extinction, not affected by the distance of the object), are shown in Fig. 22B. For the luminosity, we apply the recovered values to our M14 population (within the fiducial volume), in Fig. 22A. The resulting histogram shows three peaks between bolometric apparent magnitudes of ∼20.5, ∼21.0 and ∼21.5, corresponding to assumptions of maximum, average and minimum luminosity respectively. Realistically, the expected peak in the distribution should be located somewhere between the peaks associated with these extreme assumptions. We note that US 708 is significantly brighter (see e.g. Geier et al. 2015a) than the majority of this population, again leading credence to the idea that, if this population exists in anything approaching significant numbers, a majority of it may simply have been unobservable in the past.

Fig. 22. Expected number of objects per bolometric apparent magnitude within the fiducial volume (panel A), according to minimum, maximum and (time) averaged bolometric luminosity, and expected minimum and maximum effective temperature of He-stars in the mass range considered here (panel B). The data in (panel A) is binned along the x-axis into bins of 0.25.

9. Expected contribution to the Galactic SN Ia rate

Using our synthetic populations, we are in a position to use observations to independently calculate the birthrate of population member objects, and thus, the progenitor SN rate (which does not necessarily have to be of typical SN Ia, since the events may belong to the broader plethora of thermonuclear SNe). Specifically, as products of thermonuclear SNe, members of this putative population are, unlike the SNe themselves, or their SNRs, relatively long lived and, due to their location outside of the Galactic disc (if located within the fiducial volume), less affected by extinction, reddening and the stellar background than direct progenitor systems. Specifically, we draw attention to CD-30°11223 (Vennes et al. 2012) and the very recently reported PTF1J2238+7430 (Kupfer et al. 2022), both of which are expected to produce a runaway member of our low mass UI population with masses in the range 0.2−0.3 M_⊙, that is, proto-WDs (though in these cases the terminal accretor mass, with ∼0.9 M_⊙ is expected to be below the lowest terminal accretor mass considered here). However, we emphasise that these predictions are still, to some extent, model dependent. As such, observed products of SNe (assuming they can be claimed as such) present a model independent measure of the event rate.

A peculiarity of attempting to estimate an event rate from a population of high-velocity and hypervelocity stars is the expected inhomogeneous space density (see Fig. 6) and observational properties of the objects (sub-luminous blue, only reliably classifiable with spectra) which makes them prone to misidentification as hot WDs, O or B type stars at shorter distances or, if reddened, main sequence stars. Selection criteria focussed on high proper motion may (as discussed in the previous section) miss a large number of unbound stars and determination of radial velocities requires a larger commitment of observational resources. As the observational record may (or may not) be incomplete at the time of writing, we estimate the event rate as currently supported by observations. As, and if, additional population member are discovered, the supported event rate is likely to increase (or, less likely, to decrease, as will be explained later).

We base our calculation of the event rate on two assumptions. First, the abundance of each subpopulation (BI, BII, UI, UII) scales linearly with the ejection rate (i.e. the event rate). This assumption was tested in Sect. 6.4, and can be taken to be reasonably self-consistent. The second assumption is that the largest distance at which a population member can be detected is equal to the largest distance at which a population member has been detected. This assumption is to account for the unavailability of a reasonable estimate of the probability of detecting and correctly identifying a population member in the observational record. Further, a general assumption that individual stars, all else being equal, are more difficult to observe the more distant they are from the observer, is reasonable. We further assume the entirety of the population to be statistically distributed (i.e that the three candidate objects, as well as their distances, are no statistical aberration). As such, the farthest known population member serves as a marker of the volume used to calculate the average space density. We then use this volume to scale the space density of objects in our synthetic population to match the space density of observed objects in the same volume by linearly varying the ejection rate. The supported ejection rate can thus be written

$$\begin{array}{c}\hfill f_{\mathrm{ej},\sup }=\frac{{\rho }_{\mathrm{UI}+\mathrm{BI},\mathrm{obs}}(D_{\max })}{{\rho }_{\mathrm{pop},\mathrm{sim}}(D_{\max })}·f_{\mathrm{ej}}=\frac{N_{\mathrm{UI}+\mathrm{BI},\mathrm{obs}}(D_{\max })}{N_{\mathrm{pop},\mathrm{sim}}(D_{\max })}·f_{\mathrm{ej}},\end{array}$$(10)

where ρ_{UI + BI, obs}(D_max) = N_{UI + BI, obs}(D_max)/V(D_max) and ρ_pop, sim(D_max) = N_pop, sim(D_max)/V(D_max) are the associated number densities within the volume defined by the most distant observed object. Then N_{UI + BI, obs}(D_max) is the number of observed objects (here three: US 708, J2050, LP 398-9), N_pop, sim(D_max) is the number of expected objects, as predicted from our synthetic population, if the ejection rate is equal to f_ej = 1/300 yr⁻¹. In principle, this methodology, namely leaving the calculation of the number density dependent on the observational record (as opposed to arbitrarily choosing a volume in which to calculate the number density), has the advantage of remaining as agnostic as possible to observational uncertainties. It has, however, the disadvantage of being susceptible to further observations decreasing the supported event rate, if the next object is detected at distances significantly greater than the previous most distant object. We take this into account be treating the space density within the volume defined by the currently most distant known object, US 708, as an upper limit and the space density in the entirety of the fiducial volume (extending to D = 25 kpc) as a lower limit. In Fig. 23 we show the resulting calculated rates for the M14 population. As can be seen, the currently supported event rate, based on three tentative detections, extends from ∼3×10⁻⁷ yr⁻¹ to ∼2×10⁻⁶ yr⁻¹. We note that this value is at least two orders of magnitude lower than the predicted formation rate of likely progenitor systems for this formation channel (see Appendix A) and another order of magnitude lower than the inferred Galactic SN Ia rate (Cappellaro et al. 1997; Sullivan et al. 2006). This indicates one of the following: The probability for He-star+WD binaries to form a SN, leading to the ejection of a high- or hypervelocity star is very small indeed. This would indicate that other outcomes, such as double degenerate mergers, with the He-star evolving into a WD prior to RLOF, are favoured for this type of system. In this case, comparison of the supported rate as calculated here, with the formation rate of other possible evolutionary products, should be performed. Another possibility is that observational uncertainties have precluded detection of further population members, in which case future observations may change this picture significantly. We note that, as seen in Fig. 23, if all possible progenitor systems result in the ejection of a high- or hypervelocity He-star, we would expect of the order of 200 objects within the volume defined by the current position of US 708 (implying a detection bias of about 1/60) and of the order of 2000 objects within the 25 kpc fiducial volume. The entirety of all Galactic SNe Ia resulting in the ejection of a high- or hypervelocity He-star would result in a population exceeding 2000 objects within the volume defined by the current position of US 708 (implying a detection bias exceeding 1/600). We reasonably speculate that, even including future detections, the observational record is highly unlikely to support either the progenitor formation rate or the inferred Galactic SN Ia rate.

Fig. 23. Calculated event rates for the production of population members. Dashed lines indicate the distance-dependent event rate based on the three observed candidates. Solid lines indicate event rates under the assumption of the indicated number of detections. Distances to the observed objects (including error bars) and associated event rates are indicated by filled and hollow circles. The error bars for LP 398-9 are smaller than the radius of the circles. Upper and lower limits for the observationally supported event rate, predicted formation rate of putative progenitor systems and the inferred Galactic SN Ia rate are as indicated. For details on calculation see text.

Another possibility arises from the result that the inferred ejection rate is lower by two orders of magnitude, than the calculated formation rate of putative progenitor binaries. Namely, only a small fraction of He-star+WD binaries produce a runaway He-star. Alternatively, the production rate of putative progenitor binaries may be overestimated. However, even if the ejection rate is indeed as low as we calculate here, and if the three putative objects are no statistical aberration, then we would expect at least of the order of 25 objects (22 of them currently undetected) within the fiducial volume.

10. Conclusions

Using ejection velocity spectra for high- and hypervelocity runaway stars as presented by Neunteufel (2020), we construct a synthetic population of SN products generated in the US 708-channel. We calculate 1×10⁶ ejection trajectories each for one given terminal accretor mass, assuming ejection from the Galactic plane (i.e. z_Gal = 0) with the ejection locations determined by the Galactic density distribution. Analysing the thus generated population, we constrain the defining characteristics of members of this putative population, as well as its structure. We find that this population is composed of stars both bound and unbound to the Galactic potential. Runaways ejected while identifiable as He-stars (He-sdOs) or proto-WDs will eventually evolve into unbound or high-velocity halo WDs. We distinguish four distinct subpopulations, separated according to two independent attributes: bound and unbound (B and U), pre-WD and WD (I and II), where the UI+BI subpopulation is split along low mass (< 0.3 M_⊙) and high mass (> 0.3 M_⊙) lines, where the low mass UI+BI subpopulation will present as a proto-WD, while the high mass UI+BI subpopulation will present as a core He burning star. Thus, if the He-star+WD channel is responsible for any appreciable fraction of observable SNe, we predict the existence of a population composed of high- and hypervelocity He-stars and WDs.

In order to take account of observational constraints affecting hot subdwarf stars, we define a fiducial volume (Fig. 9), excluding most of of the Galactic disc and bulge regions. Analysing our synthetic population at the end of a 300 Myr time frame, we calculate the relative numbers of the expected subpopulations. We find that, within this fiducial volume, the profile of the BI/UI-ratio is directly related to the terminal mass of the accretor. Assuming that the channel responsible for the production of US 708 favours a single terminal accretor mass but a different terminal donor masses, this represents a directly accessible observable. We further investigate the astrometric properties of our synthetic populations at the end of the 300 Myr simulation time frame. We find that three currently observed objects, US 708/SDSS J093320.86+441705.4, as a member of the UI subpopulation, Gaia DR2 1798008584396457088 (LP 398-9), also as a member of the UI subpopulation, and SDSS J205030.39-061957.8 (J2050), as a member of the BI subpopulation, match the expected properties of the population, the former being the prototype of the population, the latter two matching sufficiently well to be designated as candidates. Of these, both US 708 and LP 398-9 are generally compatible with ejection from either Chandrasekhar and super-Chandrasekhar mass events.

We further find that the relative numbers of the population are, to good approximation, independent of the assumed ejection rate. We use this property of our synthetic population to constrain the currently observationally supported ejection rate at ∼3×10⁻⁷ yr⁻¹ to ∼2×10⁻⁶ yr⁻¹, which is at least two orders of magnitude lower than the (theoretically predicted) formation rate of the putative progenitor binary and three orders of magnitude smaller than the inferred Galactic SN Ia rate. We conclude from this that the current observational record strongly disfavours the SN channel responsible for the creation of US 708 being a major contributor to SNe Ia in general.

We further conclude that the current observational records disfavours He-star+WD binaries generally producing a member of this population. However, we emphasise that this result cannot be called conclusive in the general sense, as our models further indicate that all of the currently known objects seem to be extreme members of their own population, US 708 and J2050 exhibiting extreme radial velocities, LP 398-9 exhibiting extreme proper motion. We conclude from this that, while the extreme nature of these objects may have been a contributing factor in their original discovery over more ordinary population members, a larger population of more ordinary population members may as of yet be undetected.

If the currently low ejection rate is taken at face value, a significant number of additional, possibly detectable, population members likely exist within the fiducial volume and outside, with a significantly smaller likelihood of detection. However, should the He-star+WD binary channel for thermonuclear SNe be responsible for any appreciable fraction of the observed thermonuclear SN rate, a much larger population of runaway He-stars such as He-sdBOs (e.g. US 708 and J2050) and proto-WDs(e.g. LP 398-9), should exist, thus providing an immediately accessible observable of the associated rates. The study of high and hypervelocity stars ejected by channels not related to interaction with a SMBH should be considered a promising avenue to probe a number of unresolved astrophysical phenomena.

Acknowledgments

We thank the anonymous referee for insightful comments that improved the legibility of the manuscript. PN would like to thank Vedant Chandra for drawing attention an error in nomenclature. MK is partly supported by Grant No 11521303, 11733008, 12090040, and 12090043 of the Natural Science Foundation of China. SJ acknowledges funding from the Netherlands Organisation for Scientific Research (NWO), as part of the Vidi research program BinWaves (project number 639.042.728, PI: de Mink). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy, (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

References

Appendix A: Population synthesis

We used the fast binary population synthesis code COMBINE (Kruckow et al. 2018; Kruckow 2018) in order to calculate the formation rate of putative progenitors of HVS as under consideration in this study. The same settings as given in Table 2 of Kruckow et al. (2018) are used, excepting a constraint of the initial stellar mass range to [1 : 20]  M_⊙ for the first run (R1). To form a binary consisting of a WD and a He-star, the binary had to undergo at least one, but usually several, episodes of mass transfer. In a second run (R2) we restrict the initial mass range further to [1 : 10]  M_⊙ and assume stable mass transfer (α_RLO = 0.01, and β_min = 0.0) to be more conservative, which is more appropriate for low mass stars. Finally, we allow more efficient common envelope evolution (α_CE = 1.0 and α_TH = 1.0) in a third run (R3).

For each run, we calculate the expected formation rates in a Milky Way (MW)-like galaxy. We follow the same procedure as in Kruckow et al. (2018) by assuming a constant star formation rate of one binary per year with a primary star mass > 0.8  M_⊙ (Hurley et al. 2002). Systems of interest are selected by conditions on the mass of the He-star, M_He, the total binary mass, M_WD + M_He, and the merger timescale by pure gravitational wave radiation, t_GW. The last criterion, C8, is appropriate for the study presented in this paper. Table A.1 summarises the resulting formation rates of binaries consisting of a WD with a He-star companion. The range of 1.74 to 5.79 × 10⁻⁴ yr⁻¹ shown in Fig. 23 arises from the lower and upper limits for the full set of population synthesis runs.

All Tables

All Figures

	Fig. 1. Assumed lifetimes τ of the ejected He-stars. Stars in the shaded region are below the helium-burning limit and are deemed to thermally collapse to form WDs within 10 Myr.
In the text

	Fig. 2. Illustration of 75 tracks, out of 1×106, chosen at random, from the MWD,f = 1.4 M⊙-population. Coordinates are Galactic Cartesian. Origins, all set at z = 0, are indicated by red circles. Tracks dropping below the Galactic plane are indicated by dashed lines.
In the text

Fig. 3. Trajectories of all stars crossing into Earth local space from the outside. Unbound trajectories are indicated as red lines, bound trajectories as blue lines. Panels A and B: galactic Cartesian coordinates in the x–y-plane, with panel A showing the full trajectories and panel B the respective point of origin. Panel C: full trajectories in the x − z-plane. The positions of Earth and the Galactic centre are as indicated. Note that, as seen from above the Galactic plane (with Earth located above the plane), the Galaxy rotates in a clockwise fashion.

In the text

Fig. 4. Chirality of the ejected stars crossing into Earth local space in the M14 data set with respect to Galactic rotation. On the x-axis (location), prograde and retrograde refer to the ejection location of the star being located either prograde or retrograde in the Galactic disc with respect to Earth. On the y-axis (direction), prograde and retrograde refer to the initial velocity vector of the ejected star being oriented either prograde or retrograde with respect to Galactic rotation. The position of Earth is as indicated. Blue dots imply bound stars, red unbound stars.

In the text

	Fig. 5. Same as Fig. 3 for the US 708 local space. Here the inferred current positions of US 708 and its ejection location are indicated as well.
In the text

Fig. 6. Expected number density distribution per kpc−3 in cylindrical galactocentric coordinates for different terminal accretor masses with MWD,f = 1.0 M⊙, MWD,f = 1.2 M⊙, MWD,f = 1.4 M⊙ and MWD,f = 1.5 M⊙ depicted in panels A–D respectively. Here, the colour scale indicates the calculated number density. The current positions of Earth and US 708 are indicated. Black lines indicate surfaces of constant baryonic matter density, logarithmic in base 10, as proposed by Irrgang et al. (2013), model 1. White areas do not contain ejected stars in this realisation of our model.

In the text

	Fig. 7. Comparative number densities in the vicinity of Earth and US 708 as described in Table 4. Error bars indicate one σ.
In the text

	Fig. 8. Spectral distribution of subpopulation membership over terminal donor mass for the M10, M12, M14 and M15 populations for the total set (not merely the fiducial volume). Total numbers per mass at measured against the left hand y-axis. The black line indicates the BI/UI and is measured against the right hand y-axis.
In the text

Fig. 9. Illustration of the fiducial volume chosen for the calculation of the expected detection and event rates. Only objects located in the red shaded area are taken into account. Specifically, all objects in a distance of D < 2 kpc are assumed to be visible. Any objects at a distance D > 25 kpc are assumed to be undetectable. Objects located in the disc are discounted by the assumption that objects with a Galactic latitude α < b < −α are invisible, with α = 10°. Earth’s position is as indicated. Illustration not to scale.

In the text

Fig. 10. Cumulative distance distribution of runaway stars normalised to the total number of runaways in a non-degenerate state, Ntot, n − d for a ejection frequency of fej  =  1/300 yr−1. The blue and orange curves add up to unity, with blue indicating the relative abundance of the UI subpopulation and orange represents the BI subpopulation. Red indicates the BII subpopulation normalised to the UI+UII subpopulation. Dark orange and dark red indicate the admixture of objects in the mass range 0.7−0.8 M⊙. Panels A–D: M10, M12, M14 and M15 populations respectively.

In the text

	Fig. 11. Relative occurrence of individual subpopulations in all synthetic populations in the fiducial volume depending on different terminal accretor masses. BII/BI(τ/2) was calculated under the assumption of a 50% reduced nuclear burning lifetime of each ejected star.
In the text

	Fig. 12. Same as Fig. 10 with fej = 1/3 × 10−4 yr−1.
In the text

	Fig. 13. Like Fig. 8 but restricted to the fiducial volume. Panel A: the BI/UI ratio was cropped in order to preserve legibility. The missing value is NBI/NUI = 23 at M⋆, run = 1.00 M⊙.
In the text

	Fig. 14. Direct comparison of the BI/UI-ratio as presented in Fig. 13.
In the text

	Fig. 15. Like Fig. 13 but with assumed lifetime τ = τ/2. Panel A: the BI/UI ratio was cropped in order to preserve legibility. The missing value is NBI/NUI = 55 at M⋆, run = 1.00 M⊙.
In the text

	Fig. 16. Direct comparison of the BI/UI-ratio as presented in Fig. 15. Grey lines indicate the data presented in Fig. 14.
In the text

	Fig. 17. Like Fig. 13C, but assuming a weighting function according to Eq. (9).
In the text

Fig. 18. Spectral distribution of predicted radial velocities of the predicted UI+BI subpopulation for different terminal accretor masses within the fiducial volume. Lines of different colours indicate the sum over all runaways larger than the mass indicated in the key. Terminal accretor masses for each panel are Mf, WD = 1.0 M⊙ for panel A, Mf, WD = 1.2 M⊙ for panel B, Mf, WD = 1.4 M⊙ for panel C and Mf, WD = 1.5 M⊙ for panel D.

In the text

	Fig. 19. Like Fig. 18, but for the UII+BII subpopulation. Each line indicates the sum over all masses higher than indicated in the key; The secondary peak at about 750 km s−1 indicates the peak of the low mass UII population. The lifetimes of stars in the range 0.3−0.45 M⊙ are too long to contribute to this plot.
In the text

Fig. 20. Spectral diagram of the total proper motion of the predicted population for the M14 population. Lines of different colours indicate different distance-limited sub-samples. Dashed lines indicate where the distance-limited sub-sample of the same colour is the same as the next smaller limited sub-sample. Grey areas indicate the proper motion of the three population candidates as labelled. The PMs of US 708 and J2050 relatively similar, so their corresponding lines were unified for legibility. Proper motion errors are of the order of ∼ ± 0.2 mas yr−1 for US 708 and J2050 and ∼ ± 0.1 mas yr−1 for LP 398-9 (Gaia Collaboration 2021).

In the text

Fig. 21. Galactocentric space velocity of I-type subpopulation in the M14 (MWD, f = 1.4 M⊙) population in the fiducial volume over distance D from Earth. Colours indicate the mass of the ejected object. Panel A: entire population. Panels B and C: only contain the UI and BI subpopulation respectively. Panel D: the same as panel A, except assuming τlifetime = 0.5 τlifetime, bm. Grey indicates (with no differentiation according to mass) the II-type subpopulation in panels A and D. Positions of the observed candidates in the parameter space are as indicated. The dashed line indicates D = 25 kpc (i.e. the parameter space left of the dashed line is equivalent to the fiducial volume as defined in Fig. 9).

In the text

	Fig. 22. Expected number of objects per bolometric apparent magnitude within the fiducial volume (panel A), according to minimum, maximum and (time) averaged bolometric luminosity, and expected minimum and maximum effective temperature of He-stars in the mass range considered here (panel B). The data in (panel A) is binned along the x-axis into bins of 0.25.
In the text

Fig. 23. Calculated event rates for the production of population members. Dashed lines indicate the distance-dependent event rate based on the three observed candidates. Solid lines indicate event rates under the assumption of the indicated number of detections. Distances to the observed objects (including error bars) and associated event rates are indicated by filled and hollow circles. The error bars for LP 398-9 are smaller than the radius of the circles. Upper and lower limits for the observationally supported event rate, predicted formation rate of putative progenitor systems and the inferred Galactic SN Ia rate are as indicated. For details on calculation see text.

In the text