Predicted spatial and velocity distributions of ejected companion stars of helium accretion-induced thermonuclear supernovae

Thermonuclear supernovae (SNe), a subset of which are the highly important SNe Type\,Ia, remain one of the more poorly understood phenomena known to modern astrophysics. In recent years, the single degenerate helium (He) donor channel, where a white dwarf star (WD) accretes He-rich matter from a hydrogen-depleted companion, has emerged as a promising candidate progenitor scenario for these events. An unresolved question in this scenario is the fate of the companion star, which would be evident as a runaway hot subdwarf (He sdO/B) in the aftermath of the SN event. Previous studies have shown that the kinematic properties of an ejected companion provide an opportunity to closer examination of the properties of an SN progenitor system. However, with the number of observed objects not matching predictions by theory, the viability of this mechanism is called into question. In this study, we first synthesize a population of companion stars ejected by the aforementioned mechanism, taking into account predicted ejection velocities, inferred population density in the Galactic (Gal.) mass distribution and subsequent kinematics in the Gal. potential. We then discuss the astrometric properties of this population. We present $10^{6}$ individual ejection trajectories, numerically computed with a newly developed, lightweight simulation framework. A peak in the density distribution for close objects is expected in the direction of the Gal. center. If the entire considered mass range is realized, the radial velocity distribution should show a peak at 500\kms. If only close US\,708 analogues are considered, there should be a peak at ($\sim750-850$)\kms. We show that the puzzling lack of confirmed surviving companion stars of thermonuclear SNe, though possibly an observation-related selection effect, may indicate a selection against high mass donors in the SD He donor channel. (-abridged-)


Introduction
The study and theoretical characterization of thermonuclear supernovae (SNe) is a long-standing issue in modern astrophysics. Chief among these are SNe of Type Ia, which, distinguished by their relatively slight variation in brightness between different events, have been used successfully to measure the rate and acceleration of cosmic expansion (Nobel Prize 2011, see Perlmutter et al. 1999). Observationally, this class of stellar explosions is characterized by lack of evidence for the presence of hydrogen and is associated with a number of subtypes. These subtypes include SNe of Type Iax (Foley et al. 2013) and, possibly, certain calcium-enhanced Type Ib (Waldman et al. 2011;Neunteufel et al. 2017). Thermonuclear SNe derive their name from their assumed origin as thermonuclear explosions of white dwarf stars (WDs), the degenerate remnants of low to lowintermediate mass ( 8 M ⊙ ) main sequence (MS) stars. Here, detonation is triggered by mass transfer from a companion star of some description. This explosion mechanism is distinguished, on the theoretical side, from core collapse SNe, generally triggered by the gravitational collapse of the core of higher-mass (8 M ⊙ M 40 M ⊙ ) MS stars (see Sukhbold et al. 2016, for a review) and, at the highest masses, pair instability SNe become a possibility (see Heger & Woosley 2002, for reviews). Note that PISN are also thermonuclear in nature, but without involvement of a WD.
Thermonuclear SNe (of the type under consideration here) are generally categorized, depending on the structure of their progenitor systems, into two distinct channels: Double degenerate (DD), where the companion is another WD, and single degenerate (SD), where the companion is a non-degenerate star, such as a MS star (see Hillebrandt & Niemeyer 2000;Ruiter 2020, for reviews).
An open question in studies of the SD channel is that of the nature of the non-degenerate companion, with stars on the MS, red giant branch (RGB), hydrogen depleted stars (hereafter: He-stars) either during or after their core helium (He) burning phase having been proposed and studied by different authors (e.g. Nomoto 1982; Woosley & Kasen 2011;Brooks et al. 2017).
However, regardless its structure, the donor star is expected to survive the accretor's detonation and, hence, the SN, in all viable SD scenarios (see, e.g. Liu et al. 2015;Liu & Zeng 2020). If, in the aftermath of a thermonuclear SN, the accretor is com-A&A proofs: manuscript no. HVS-pop pletely disrupted, then the former donor star is, according to conservation of orbital momentum, flung away from the location of the progenitor system at a velocity, due to ejecta interaction, slightly higher than its pre-detonation orbital velocity, as dictated by the involved masses and the binary orbital separation (Bauer et al. 2019). Discarding certain super-Chandrasekharmass mechanisms (see, e.g. Di Stefano et al. 2011), which lie outside the parameter space considered in this study, detonation is presumably triggered by mass transfer from the nondegenerate donor star, the donor star can be assumed to be filling its Roche lobe at the point of detonation. The final attainable orbital separation, and thus runaway velocity, is thus dictated by the radius of the donor star. Consequently, the lowest runaway velocities are expected in progenitor binaries containing a low mass RGB donor, followed by He-giants and typical MS stars. The highest runaway velocities in the SD channel are expected to require a low-mass, hydrogen depleted star (Justham et al. 2009;Neunteufel 2020).
For most reasonable terminal mass combinations, the ejected donor star in this scenario is expected to move at sufficiently high velocities post-ejection to become unbound from the Galaxy. Such unbound stars are, in addition to their other attributes, commonly called hypervelocity stars (HVS). While many observed HVS, particularly the more massive specimens (Hills 1988;Brown 2015, the latter being a recent review) are thought to originate from encounters of binary star systems with the supermassive black hole at the center of the Galaxy, with one object (Koposov et al. 2020) having been traced back to this origin with high probability. HVS originating from the Galactic (Gal.) center have also attracted attention as tracers of the Gal. mass distribution (Yu & Madau 2007). A substantial body of work (e.g. Brown et al. 2009;Lu et al. 2010;Rossi et al. 2014) has considered the kinematic properties of these objects. Possible origins as former companion stars have also been considered (Evans et al. 2020) and the search for HVS of all masses is ongoing (Raddi et al. 2020;. However, at the low-mass end of the mass distribution, a number of objects have been detected (Geier et al. 2015;Vennes et al. 2017;Shen et al. 2018), likely ejected in thermonuclear SN events. Of these, the hot subdwarf US 708 is of particular interest, as its current structure, composition and velocity suggest it as a likely product of the He donor SD channel (Geier et al. 2015;Neunteufel 2020, hereafter G15 and N20 respectively).
Unlike surviving companions resulting from different candidate progenitor channels, US 708 analogues (i.e. hypervelocity hot subdwarfs) are, again due to their composition and velocity, relatively clearly identifiable, if not easily detectable. However, with direct observation of clear thermonuclear SN progenitor systems not forthcoming, analysis of their surviving ejected companions would offer an avenue of investigation into the terminal state of thermonuclear SN progenitor systems not available otherwise. An open question in this regard is whether the serendipitous detection of US 708, the only known US 708 analogue at the time of writing, indeed suggests a, as of yet undetected, population of these objects. Belonging to the class of He-enriched hot subdwarf O stars, US 708 is currently the only known unbound hypervelocity He-star. To date about 6000 hot subdwarfs are known (Geier 2020;Heber 2016), which are found in all Galactic populations (Luo et al. 2020). Several fast, yet still bound hot subdwarfs have been discovered (Tillich et al. 2011;Németh et al. 2016;Ziegerer et al. 2017) and the halo population contains many more sdO/Bs with quite extreme orbits (Luo et al. 2020). It is therefore not straightforward to distinguish between bound runaways and halo stars based on their kinematics. This paper is aimed at predicting the most likely kinematic properties of such a hypothetical population of US 708 analogues, taking into account, as much as possible, uncertainties in the areas of explosion mechanism, terminal component mass and predicted scatter in runaway velocities. Movement in the Gal. potential and the Gal. density distribution are fully considered.

Physical assumptions
As stated above, the SN channel under consideration here is the SD He-donor channel. More specifically, we consider the case of the mass donor being a low-mass hydrogen depleted, core-He burning star (i.e. a He sdO/B). This channel was studied in detail by N20, who showed that the ejection velocities of runaways produced with this mechanism are largely dictated by the terminal donor and accretor masses, not the history of the binary following the emergence of its terminal structure as a system consisting of a He subdwarf and a WD, nor the amount of transferred material. The ejection locations in this study are assumed to lie in the Gal. disk and outside of clusters and close associations, but we do not expect our results to change significantly if this assumption is relaxed (See detailed discussion in App. B).
This study relies on a combination of the results of detailed stellar evolution calculations with stellar kinematics (as detailed below), which does allow a large number of experiments to be performed while being computationally relatively inexpensive.

Initial conditions of the run
The initial conditions for our population consist of the ejection velocity for each ejected runaway subdwarf, which is a function of its mass, the ejection location, i.e. the location of the progenitor binary within the Gal. mass distribution, and the ejection direction, i.e. the direction, with respect to Gal. rest frame coordinates, in which the runaway is ejected. Each of these parameters is generated randomly (see App. A.2). We briefly note that velocities, as a function of mass, are chosen based on the velocity spectra presented by N20 (see App. A.3 for details), all ejections are assumed to originate in the Gal. plane (i.e z Gal = 0) and the initial radial distance from the Gal. center is weighted (see App. A.4 for details) according to the density distribution obtained by Irrgang et al. (2013).

Kinematic simulation
Based on our randomly generated set of initial conditions, we use the SHyRT framework (See App. A.1) to calculate the resulting runaway trajectories. Specifically, using the given model Gal. potential, we calculate the average orbital velocity at the point of ejection and, taking into account the ejection direction and ejection velocity with respect to the center of mass of the progenitor binary, we calculate the Gal. rest frame velocity of the ejected star. We then allow each ejected star to move in the Gal. potential (see App. A.5 for details) for a span of 300 Myr, which is somewhat higher than the maximum nuclear timescale of a low mass (∼ 0.45 M ⊙ ) He-star (see App. A.6).

Synthesis of the runaway population
The resulting trajectory data is then post-processed to generate the expected runaway population. This is accomplished by assigning each trajectory a timestamp corresponding to time intervals mimicking the inferred Gal. SN Ia rate. The assumed inter-val is chosen to be 3 · 10 −3 yr −1 (Cappellaro et al. 1997). As our trajectory calculations use a basic timestep of one Myr, intermediate positions are interpolated from the data. Furthermore, trajectories are truncated according the nuclear timescales of typical He-stars of the respective masses (see App. A.6). Timedependent evolution of the population can be expected. Truncation of the He-star lifetime leads to an equilibrium between objects being removed from the sample and replenishment of the population through continuous ejection. Close US 708 analogues (with masses in the range 0.35 M ⊙ ≤ M ≤ 0.4 M ⊙ ), due to their lifetimes exceeding 300 Myr, are not affected by this prescription. These objects will be ejected with velocities exceeding 800 km s −1 and can be assumed to have moved sufficiently during that time to likely be inaccessible to observation.

Results
The reader should note that we base our analysis on two different premises: First, we assume that the population consists of runaways, in a flat distribution, of the entire range of masses (0.2 M ⊙ ≤ M ≤ 0.8 M ⊙ ) and, second, a truncated distribution with only runaways of masses 0.35 M ⊙ ≤ M ≤ 0.4 M ⊙ taken into consideration. The latter is chosen to represent close analogues of the observed star US 708, which was previously inferred to have a current mass in the range 0.3 M ⊙ ≤ M ≤ 0.4 M ⊙ (G15, N20) with a terminal accretor mass corresponding to the assumed 1.4 M ⊙ in this study. This, referring to Bauer et al. (2019) The former distribution will be referred to as the "full sample", while the latter will be referred to as the "reduced sample". The total number of objects in the reduced sample is about 15% that of the full sample. Note also that the main goal here is not to provide predictions of the total number of observable objects, as uncertainties in the production rate of the progenitor systems, as well as detonation mechanisms, inhibiting unambiguous identification of spectrally classified SNe and progenitor systems, are still too great to overcome (see, however, Wang et al. 2013, who estimate the rate of He donor induced SNe to account for about a third of the inferred Gal. SN Ia rate). Instead we present relative normalized number densities, corresponding to a probability distribution for the most likely observable characteristics.

The galactocentric distribution
We present results for the galactocentric distribution of our synthetic population in Fig. 1. As shown in Figs. 1A and B, the resulting distribution is essentially cylindrically symmetric. The apparent peak in the density distribution around the Gal. center (x Gal = 0, y Gal = 0) is a result of the, somewhat dispersed, initial distribution conforming to the Gal. density distribution. A noticeable peak in the density distribution at the location of the Gal. disc (z Gal = 0) is apparent, reflecting both the initial distribution (projected onto the z-axis) and deflection of the random ejection direction due to Gal. rotation. We conclude at this point that, neglecting instrument related selection effects, observations focused toward the Gal. center will yield a higher probability of new discoveries. Further, observations should (assuming that the entire mass range is encountered in reality), yield a higher ratio of more massive runaways close to the Gal. center and the Gal. disc, numbers of less massive (and therefore longer lived), runaways being reduced as these objects are ejected. Concurrently, more massive runaways will encounter the end of their He burning lifetime at shorter distances from their ejection locations both due to their shorter lifetimes and their slower ejection velocity.

The observational distribution
The celestial coordinate distribution in right-ascention (RA) and declination (DE) of our sample, as well as the distance, are summarized in Fig. 2. Note that here we focus on "closer" objects (Distance D < 20kpc) in both subfigures (A) and (B). We note a significant peak both in RA and DE (RA = 275 deg, DE = -27 deg). These peaks are present in both the full sample and in the reduced sample and coincide with the position of the Gal. center. Noticeably, the current position of US 708 does not correspond to the predicted positions of the peaks. The distance distributions (measured from the position of Earth) show a plateau at ∼ 20 kpc (this forms a peak when larger distances are included) in the full distribution. This is in contrast to the reduced distribution where the same peak is located at ∼ 190 kpc (See App. D). These peaks are essentially a reflection of the object number densities being integrated over the entire celestial sphere (compare Fig. A.2). However, the location of these peaks is instructive of the effect discussed in Sec. 3.1, with numbers of lower mass runaways being attenuated due to their higher runaway velocities. If the entire mass range is realized in nature, observations should encounter a relatively higher number of higher mass runaways in the Solar environment than lower mass ones. While there may be a significant observational bias against the detection of higher mass He-stars (see discussion in App. C), the only observed object being of a lower mass may be an indication that higher mass runaways are selected against. The latter conclusion is significant insofar as donor stars in the He-donor scenario for thermonuclear SNe, particularly in the double detonation scenario, are thought to possess higher masses under most conditions (M > 0.8 M ⊙ , see Neunteufel et al. 2016).
We present the distribution for parameters of motion in Fig. 3. Again we truncate the sample at a distance of D < 20 kpc. We give the parameters of motion for US 708 for reference. While there are pronounced peaks in the proper motions for the full sample, these peaks are much less pronounced in the reduced sample. In both cases these peaks are located around the zero axis, reflecting objects at greater distances dominating the population. In both cases values of more than 20 mas yr −1 can be expected, mostly representing objects being ejected from the Gal. center perpendicular to the line of sight. In both cases, the distribution being shifted towards positive velocities is a geometric effect. Noticeably, there is a pronounced peak in the radial velocities present in both the full sample (∼ 500 km −1 ) and the reduced sample (∼ 750 − 850 km −1 ). US 708 lies in the high-velocity tail of both samples. In any case, we conclude that proper motion is likely not a promising observational filter for these objects with searches for high radial velocities being more attractive.

Conclusions
Combining preexisting data obtained via detailed stellar evolution calculations, model Gal. potential and density distribution and numerical kinematic calculations, we have constructed a synthetic population of US 708 analogue HVS from a sample of 10 6 individual ejection trajectories. From this synthetic distribution, we derive predictions for the most likely observational characteristics of members of this population. Most importantly, we find that, if the entire mass range of ejected Hestar masses is encountered in nature, there should be a significant excess of higher mass runaways over lower mass ones close to Earth. This prediction is at odds with the only observed object, US 708, which is likely a lower mass star, This, in itself, contains a discrepancy at least with theoretical models on the double detonation scenario predicting higher mass donor stars. This discrepancy may imply that higher mass donor stars are selected against in the pre-SN evolution of the progenitor binary, in which case all observed runaway He donor stars can be expected to possess masses comparable to that of US 708. It is further evident that, should no new discoveries of US 708 analogues be forthcoming, the He donor SD channel for thermonuclear SNe would be significantly challenged. We further find that high velocity runaways in our sample will primarily be distinguished by high radial velocity, with the probability distribution of the entire sample peaking at ∼ 500 km s −1 , the reduced sample peaking at (∼ 750 − 850 km s −1 ), with US 708 located in the highvelocity tail of either distribution. We note that US 708 has been discovered based on a spectrum taken during the Sloan Digital Sky Survey (Hirsch et al. 2005) and most subsequent searches were based on new releases of SDSS. Hot subdwarfs are most easily detected as faint blue stars at high Galactic latitudes. According to our predictions, the currently available survey area at high Galactic latitudes in the Northern hemisphere turns out to be worst choice to search for sdO/B runaway, because it leaves out both the disk and the bulge. Southern surveys such as SDSS-V and 4MOST will likely be much more suited in this respect. We conclude that the apparent lack of US 708 analogues might to a large extend be due to strong observational selection effects. These can be overcome by conducting deep spectroscopic surveys in the South, obtain high-quality spectra of carefully selected candidates and improve the models to fit them, and also by checking the classifications and parameters of potentially misclassified bright candidates from the literature. However should future observation campaigns fail to produce further US 708 analogues the overall viability and prevalence of the SD He donor channel should be seen as seriously challenged on the theoretical side.
A&A proofs: manuscript no. HVS-pop Appendix A: Remarks on the methodology Appendix A.1: Numerical framework In this study, we use a lightweight stellar kinematics framework ("SHyRT") to compute the ejection trajectories of our sample. The "SHyRT"-framework (Simulated Hypervelocity Runaway Trajectories) is a newly developed, lightweight numerical tool designed to rapidly compute movement of ejected stars in the Gal. gravitational potential and provide synthetic astrometric observables of the calculated trajectories. At its, core, the SHyRT framework solves the basic Newtonian equations of motions in the form utilizing a fourth-order Runge-Kutta integration scheme with adaptive step size control as described by Press et al. (1992). The reliability of this tool was verified by comparison with extant tools designed with the same goal in mind, such as gal_py and is identical to that utilized by N20. Importantly, for this study, the SHyRT-framework calculates Gal. rotational velocity curves and includes them in the initial conditions for each simulated trajectory.

Appendix A.2: Initial conditions
The initial conditions for this simulation consist of seven primary parameters: Position in three dimensions, the ejection direction in three dimension and the total velocity with respect to the center of mass of the progenitor system. All of these parameters are randomized as far as possible, using the prescription described below. For the velocity, we adopt the ejection velocity spectra presented by N20, (pseudo-)randomly 1 choosing a discrete mass in the range 0.2 M ⊙ ≤ M ≤ 0.8 M ⊙ (0.05 M ⊙ intervals) for the donor/runaway, then randomly choosing a velocity appropriate for that mass according to the spread of the spectrum (see also Appendix A.3). We make a conscious choice to only consider discrete values for runaway masses for convenience in population analysis. The chosen mass range corresponds to the ejection velocity spectra available from N20, the lower limit in that study being dictated by considerations of numerical stability and the upper limit by a desire to explore the production of HVS, with 0.8 M ⊙ being close to the upper mass limit for production of HVS in the solar neighborhood. For the ejection location, we assume that all stars originate from the Gal. plane, so no variation was assumed along the Gal. z-axis (though movement is fully allowed in 3D). For the remaining two dimensions, a polar coordinate system was adopted, and an angle in the range 0 − 2 π was chosen randomly. The radius was also chosen randomly, though weighted such that the line number density of ejection locations obeys ρ N (R) ∝ ρ(R) · R, where ρ N (R) is the number density at distance from the Gal. center R and ρ(R) 2 is the number density at the same location (See App. A.4). The ejection direction is randomly chosen to conform to a flat isotropic distribution.

Appendix A.3: Ejection velocity spectra
Ejection velocities were chosen according to the ejection velocity spectra calculated by N20 ( Figure 8E) in the following way: 1 unless stated otherwise, random numbers were generated using the rd()-function, contained in the <random>-library for c++. 2 i.e. the total number of objects is N = x Gal [kpc] Sun Gal.Cen. The radial distribution of ejection location is shown in Fig. A.2. We deem the randomized distribution to adequately conform to expectations. The distribution of ejection locations in the x-y-plane is shown in Fig. A.3. We note in passing, that the most likely ejection location of the observed object US 708 lies close to the maximum in line density at ∼ 6 kpc. . If so, the ejected object will form first a proto-WD, then, after further contraction, a (extremely-) low mass WD. There is, at the time of writing insufficient data available in the literature to cover the entirety of our parameter space and all eventualities.
In order to confront these issues we adopt a the following assumptions and prescriptions: We neglect the difference in remaining He MS lifetimes of our runaways. Instead, we adopt certain remaining lifetimes according to terminal He donor mass, corresponding to the nuclear timescale of a typical He-star of the same mass. This implies that at, ejection, each runaway has the entirety of its normal He burning lifetime ahead of it (except if M ≤ 0.3 M ⊙ , as discussed below). This leads to an underestimation of each star's luminosity, as most stars of the type under consideration here are expected to brighten over time, but an overestimation of stars in their observable phase since, as mentioned above, most of the objects here will have a diminished remaining He burning lifetime at the point of ejection. The assumed lifetimes are given in Fig For helium-enriched stars such as US 708, however, the systematic uncertainties on crucial parameters such as the surface gravity are still quite high and make a precise mass determination difficult. Given these problems, the yet missing population of intermediate-mass runaways might just be hidden among the already discovered stars. And finally, intermediate-mass helium stars can be quite luminous and can therefore be detected far away. Since they are predicted to be close to the disk, the substantial reddening caused by interstellar matter in the disk will then change their apparent colours significantly and move them out of the search spaces for our survey . Furthermore, such stars can be easily misclassified as main sequence stars as in the famous case of the proposed black hole binary LB-1 (Irrgang et al. 2019).