© ESO 2022

1. Introduction

Decades ago, it was recognized that there are young stars present outside the spiral arms. These stars seemed to have left their place of birth and were moving away from the Galactic disk. Blaauw & Morgan (1954) studied the proper motions of two such stars, AE Aurigae and μ Columbae, and noticed that they were moving away from each other in opposite directions, with the direction of movement pointing away from the Orion nebula. These two stars formed a part of the larger survey of such stars by Blaauw (1961), where the name “runaways” was used if their peculiar velocities were greater than 40 km s⁻¹. Blaauw (1961) explained that such stars had been part of a binary system, and when one of the stars underwent a supernova, the other star was released with a velocity close to its original orbital velocity; we refer to this scenario as the binary supernova scenario (BSS). Many runaway stars have been discovered (e.g., Gies & Bolton 1986; Heber et al. 1997, 2008; Tenjes et al. 2001; Lynn et al. 2004; Dinçel et al. 2015; Renzo et al. 2019; Raddi et al. 2021) since then.

An alternative scenario of dynamical ejection from stellar clusters was proposed (Poveda et al. 1967; Fujii & Zwart 2011; Hoogerwerf et al. 2000). Accordingly, stars can be ejected from a stellar cluster during its initial dynamical relaxation. This scenario, known as the dynamical ejection scenario (DES) also explains the expansion of compact and loose stellar clusters (Bastian et al. 2008; Pfalzner & Kaczmarek 2013; Farias et al. 2019).

Regardless of the formation mechanism, most of the stars discovered so far are massive stars of spectral type O and B and can be divided into two types of young runaway stars. Young massive stars found in the halo which have already left the disk of the Galaxy (Silva & Napiwotzki 2011), and stars which seem to be traveling away from stellar clusters at space velocities of 30 km s⁻¹ (Gies & Bolton 1986) or more, but have yet to leave the disk. This lower velocity limit is somewhat arbitrary, but it separates the runaways from the class of walkaway stars, usually defined simply as being unbound from their parent associations (de Mink et al. 2014; Renzo et al. 2019). Walkaway stars are similar to runaways in terms of their origin, but they are more abundant due to their lower velocities (de Wit et al. 2005). Since the average velocity dispersion in such associations has been found to be between 3–8 km s⁻¹ (Mel’Nik & Dambis 2009; Mel’Nik & Dambis 2017), a strict threshold on the minimum velocity of walkaway stars has not been set. However, previous studies have used a lower value of 5 km s⁻¹ (Eldridge et al. 2011) to account for all walkaways, or 10 km s⁻¹ to take the maximum velocity dispersion in all associations into account (Schoettler et al. 2020).

Although studied for decades, the origins of most runaways are still not known. Nevertheless, interest in kinematically peculiar stars, especially with the most extreme velocities, has never been lost. Hills (1988) had theorized that stars with velocities above 1000 km s⁻¹ might also exist. A supermassive black hole at the center of the Galaxy could tidally disrupt a tight binary system and release one component at a velocity up to a few thousands km s⁻¹, which would be unbound to the Galaxy. However, it took until 2005, when hypervelocity stars that are unbound to the Galaxy were discovered (Brown et al. 2005; Hirsch et al. 2005; Edelmann et al. 2005). The Gaia mission (Gaia Collaboration 2018) provided proper motions which, for the first time, allowed us to study the candidates in the full six-dimensional phase space. Indeed, SS5-HVS1 was discovered as an A-type star traveling at a rest frame velocity of 1755 km s⁻¹ (Koposov et al. 2019). Gaia-based results for other B-type candidates indicated that many of them do not originate from the Galactic center but from the Galactic disk (Irrgang et al. 2018). Despite the high end of the velocity spectrum receiving much attention, studies have continued to show that most stars do not come from the center of the Galaxy (Irrgang et al. 2021). The suspected source of origin in this case are star clusters. After birth a few of the stars are ejected either dynamically through interaction with other stars, or released when its companion in a binary or triple system undergoes a supernova. The aim of this paper is to identify the parent clusters of known massive runaway stars (mostly of spectral type B). In Sect. 2 our calculations of trajectories in a Galactic potential is described. Section 3 describes a new fitting method. The Orion nuclear cluster is used as a testbed for the fitting method (Sect. 4). In Sect. 5 we revisit the sample of cluster runaways investigated by Hoogerwerf et al. (2001). In Sects. 7 and 8 we apply our new method to B type runaways in the Galactic halo and summarize the results and conclude in Sect. 9.

2. Trajectories

To provide evidence of origin of stars in clusters, we need to trace back stars and clusters in time. There are two main requirements to solve this problem. Firstly, the Galactic potential of the Milky Way needs to be modeled, which requires an understanding of the structure of the Milky Way. Secondly, the trajectory of a test particle (in this case a star or a cluster) under the influence of a potential needs to be calculated numerically.

Allen & Santillan (1991) proposed an analytic model for the Milky Way’s gravitational potential. Their model consists of three axis-symmetric components: A spherical central mass or bulge, a disk, and a spherical halo. The model has since been used to calculate and pinpoint the orbits of hypervelocity stars in many studies, for example by Odenkirchen & Brosche (1992). The parameters of this model have been updated by Irrgang et al. (2013) to match recent observations (e.g., Irrgang et al. 2010). We use this model (model I of Irrgang et al. 2013) for all computations of trajectories using the most recent version of the orbit integrator by Irrgang et al. (2021) as well as galpy (Bovy 2015). The second part of the analysis is to find out if two objects come from the same place at the same time, which has so far required Monte Carlo simulations. Therefore, the kinematic parameters (distances, positions, proper motions, and radial velocities) of both of the objects are first acquired.

Typically, a Monte Carlo simulation is run which randomizes the value of the parameters within their uncertainties. For each run, the trajectory of both objects in the Galactic potential is computed, and the minimum distance between the two is calculated for all possible combinations of trajectories. This process requires a computation of millions of trajectories and the acquired histograms provide insights into whether the two stars (or a star and a cluster) were close enough to be considered as having a similar origin or not. In this paper, we shall employ a new method of finding sites of origin for runaway stars. Our method is based on a minimization procedure, and we shall use it on candidate lists of stars provided in Sects. 5–8.

3. Fitting procedure

Hoogerwerf et al. (2001) provide mathematical evidence for why it is not possible for the distribution of the closest distance of encounter between two objects to peak at zero, even if the two objects come from the same space-time point. The reason is that in three-dimensional space the errors on the observed parameters shift the observed mean of the closest possible distance. For two objects we have at least eight parameters (the two parallaxes, the four proper motion components, and the two radial velocities) with nonzero standard deviations. To claim that two objects come from the same place at some time, they relied primarily on simulations of trajectories of the two stars in a Galactic potential. These simulations have since been used by other authors (Irrgang et al. 2010; Neuhäuser et al. 2020a; Bobylev & Bajkova 2021) to search for the origins of runaways and hypervelocity stars. However, these simulations have two shortcomings. Firstly, they are time consuming. For example, to model only a few thousand combinations of trajectories imposing an encounter distance of less than 10 pc, at least a million combinations of trajectories need to be calculated. For instance, for ζ Oph and PSR J1932+1059, only 4214 trajectories out of three million trajectories simulated by Hoogerwerf et al. (2001) came closer than ten pc. Secondly, due to time and memory constraints, simulations are highly inefficient for shortlisting candidate lists of potential runaway stars. This is especially pertinent for candidates which are more than a kpc away, and might require computations of trajectories up to hundred Myr.

In the following we introduce a way to more efficiently calculate whether two or more objects can be said to come from the same place. Our method is based on a fitting procedure which fits the trajectories of the objects by minimizing the distance between them using a χ² method. This technique is computationally much less costly than Monte Carlo simulations.

The fitting procedure allows us to find the minimum distance possible between the trajectories of two or more objects. The trajectories are calculated using a gravitational potential, as defined in Sect. 2. The first trajectories use the mean values of the kinematic parameters. We then calculate the distances between all points in these set of trajectories, which allows us to find the minimum distance for the set of trajectories. This distance is then required to be minimized in order to find trajectories which have the least minimum possible distance. The complete algorithm then followed the following steps:

1. The trajectories for two or more objects were computed using the mean values of kinematic parameters (for n objects these are 6*n mean values). 2. The distance between the objects was calculated for all time steps till the maximum allowed time. The minimum of this distance provides the first guess for the minimization procedure. 3. This distance was kept as the function to be minimized. This was followed by minimization of the distance function, wherein the kinematic parameters were allowed to vary within their uncertainties (for n objects these were 4*n parameters, since the initial right ascension and declination were kept constant). 4. The minimization was stopped when a certain lowest distance was reached. Ideally, the distance for two point like sources may be close to zero, although the required threshold was kept higher.

A drawback of the χ² distance minimization approach is that the trajectories may be statistically insignificant. For instance, it is possible to find trajectories which have parameter values 4–5σ away from their means, simply due to the fact that these trajectories might provide the best fits. Therefore, to remedy this, we created another measure for our minimization procedure which took into account the significance of our trajectories by weighing the significance of the final trajectory in our distance function through the use of the statistical p-value.

We calculate the average distance between ‘N’ objects at any time as a function D_N. This distance (given here in kiloparsec) is calculated using the Galactic potential model 1 of Irrgang et al. (2013). For more than 2 objects, we take the average distance between them. To add a measure of significance, we define the parameter D_p as:

$$\begin{array}{c}\hfill D_p=\mathrm{wf}\ast (1-p).\end{array}$$(1)

The multiplicative factor wf is a weighting factor which allows us to control how much weight we want to allow this function to have. In our algorithm we use the p-value, given as p here, which we calculate using the χ² method. Since we have many parameters, the final p-value is calculated using all the parameters and serves as a metric for how close the trajectories are to the mean trajectories.

Finally, we minimize a function which takes into account both D_N and D_p, defined as:

$$\begin{array}{c}\hfill D_M:\sqrt{D_N}+D_p.\end{array}$$(2)

The final value of D_M is therefore highly dependent on the choice of weighting factor. We show here that subject to the condition that D_N be less than a predefined chosen value δ, it is possible to choose a wf such that both D_N and D_p are almost equally weighted. The chosen value δ then defines a ‘close encounter’. For our fitting procedure throughout this paper we choose δ as 1 pc¹. Therefore, considering this value to be the maximum possible minimized value, we choose to have D_N = 0.1 pc to validate that two or more objects come from the same place. We assume that the p-value should be roughly 0.9. For reference, a value p > 0.43 for 8-parameters corresponds to all 8 parameter values being within 1σ of their means. A value of 1.0 may be chosen, but is practically very difficult to achieve when performing any minimization. Therefore, for equal weighting we have from Eqs. (2) and (1),

$$\begin{array}{c}\hfill \sqrt{{10}^{-4}}={10}^{-1}\ast \mathrm{w}.\mathrm{f}\Rightarrow \mathrm{w}.\mathrm{f}={10}^{-1}.\end{array}$$(3)

Hence, our procedure minimizes the quantity:

$$\begin{array}{c}\hfill D_M=\sqrt{D_N}+0.1\ast (1-p).\end{array}$$(4)

This implies that step 4 in our procedure was changed to: Stop the minimization when the best fit is found for a distance lower than the threshold distance and a p-value higher than the chosen significance level.

The adopted null hypothesis states that the two or more objects whose trajectories we are optimizing have an optimal trajectory which originates from the same space-time point. The significance level, a p-value below which would allow us to reject the null hypothesis, should therefore be decided before hand. This level, also referred to as α, depends on the number of parameters we fit for our trajectories, as well as the method used to compute it. Therefore, common α values are not conceptually the correct ones for our work, since they are chosen based on 1 parameter Gaussian distributions. As can be seen in Table 1, the 8 parameter minimization can lead to very low p-values when compared to the 1 parameter solution. For our present paper, we chose the 2σ rejection limit for 8 parameters. Actually, we aim at the failure of rejecting the null hypothesis and, therefore, we regard any p-value above the 1σ rejection limit as support to the null hypothesis. This implies a stringent condition that there is a common origin for two objects where all parameters of the objects are within 1σ of their observed means.

3.1. Extent of clusters

The minimum spatial distance δ, considering the extent of open clusters on the sky, may be considered acceptable up to 10 pc for close encounters. However, we are more conservative by requiring that the minimum distance be less than 1 pc from the cluster center. Most of the runaway ejections are thought to happen here due to the density and compactness of the cluster cores. Cantat-Gaudin & Anders (2020) provide a list of clusters with parallaxes and angular extent of the clusters. We used these to compute the radius of the clusters within which 50% of the cluster members are located (see Table 2 and Fig. 1). The mean of the histogram is ≈3.7 pc with a standard deviation of 2.4 pc. This is slightly higher than our chosen significance level and therefore, in case the minimization fails for a particular star/cluster or association, the cases are inspected individually allowing for somewhat larger distances for close encounters. For instance, the stars HIP 22061 and HIP 29678 both come close to the cluster Collinder 69, that is 7 pc from the center of the cluster, which is identical with the extent of this cluster. Hence, the two stars may have been ejected from the outskirts of the cluster rather than its central region. The rigid criterion also allows us to include small clusters.

Fig. 1. Cluster extent.

Table 2 shows the radius for some of the new clusters we have studied in this paper.

3.2. Age of clusters and runaway stars

The high accuracy of Gaia astrometry has allowed the CMDs of open clusters to be cleaned by improving cluster membership from precise common proper motions. Hence, the cluster proper motions are based on large number of cluster members (see Table A.2). Based on such cleaned CMDs ages were recently derived by isochrone fitting for many clusters (e.g., Bossini et al. 2019; Liu & Pang 2019). In many cases cluster ages are given in the literature with statistical uncertainties only, based on the favored isochrone models. The role of different stellar models has been investigated in a few cases (e.g., for the well studied cluster Alessi 13, Galli et al. 2021). While common origin of two or more objects can be constrained spatio-temporally by using the distance between them at a certain time, the characterization of the method of origin into BSS or DES requires further study. The DES method happens during the formative years of the clusters (Hoogerwerf et al. 2001; Fujii & Zwart 2011). Therefore, for early type massive stars the cluster age and flight time should be comparable. On the other hand the BSS requires the evolution time of the primary star. This implies that the cluster should not be so old that it precedes the formation of the O- or B-type stars that were the primary, and the difference between the age of the cluster and the time of flight of a BSS runaway should be bounded. For this reason, a cluster may theoretically have a maximum of age of ≈150 Myr (Georgy et al. 2013), which is roughly the age of a 4 M_⊙ B-type star. This age can be slightly greater, since BSS candidates can appear younger due to the rejuvenation by the primary. They are then known as blue stragglers, which have been predicted to exist as a result of collisions in open clusters (Portegies Zwart et al. 1999), have been expected and found in ejected populations of runaway stars (Portegies Zwart 2000; Hoogerwerf et al. 2001), and have even been suggested to explain the place of origin of hypervelocity stars (Perets 2009; Brown et al. 2010). Theoretically, considering a minimum mass of nine solar masses required for a supernova, the approximate lifetime of 30 Myr (Georgy et al. 2013) should be considered for a BSS progenitor, and therefore the difference in the cluster age and the flight time of the runaway star should not exceed 30 Myr. Due to the existence of blue stragglers and the fact that clusters ages might be inaccurate, this difference may appear greater in practice.

In the following we first look for clusters with a spatial coincidence with runaway stars in the past. We then use cluster ages to confirm or rule out common origins.

3.3. A comparison of measurements

To check whether the accuracy and precision of the new measurements is comparable to our distance thresholds, we ran 100 000 simulations of orbits for AE Aurigae and μ Col, using HIPPARCOS second reduction and Gaia EDR3 data. The results shown in Fig. 2 for EDR3 show most distances close to 10 pc (lower panel), which is more than an order of magnitude improvement compared with results from HIPPARCOS data (upper panel and is mainly a consequence of the factor 10 improvement of the proper motions (see Fig. 3). Owing to this, we find that our χ² fitting procedure, which is designed to target significant orbits with closest distances below 1 pc, may indeed be used with the present day measurements.

Fig. 2. Histograms for back-traced orbits of AE Aurigae and μ Col from Monte Carlo simulations. with HIPPARCOS second reduction (top) and Gaia EDR3 data (bottom). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative velocity of the stars with respect to each other.

3.4. AE Aur and μ Col

The first tests done to validate our method were for the two well studied stars AE Aurigae and μ Columbae, which have been shown to have originated from the same space ≈2.5 Myr ago. We relied on the same HIPPARCOS astrometry as Hoogerwerf et al. (2001) to carry out a first test. Using the minimization method above and allowing for a variation of each parameter up to 2σ uncertainty, allowed us to find the best fit times and the corresponding p-values. We found a p-value of 0.72 using HIPPARCOS astrometry, suggesting a significant encounter between the two stars 2 Myr ago. This p-value was even higher when the parameters were allowed to vary only till 1σ uncertainty, which is a much more stringent condition in this scenario. We found a p-value of 0.87, which implies that there exist trajectories for which AE Aurigae and μ Col originated at the same point, with all kinematic parameters close to their observed means.

As a test of our metric, we also ran the fitting but without the p-value term (wf was kept 0). For the 2σ case, the a posteriori calculated p-value of the minimized trajectory was 0.002, a drastic reduction from the value of 0.72 before. This allows us to conclude that the addition of the p-value term into the distance metric has two important consequences. Firstly, the explicit p-value term allows to find more statistically significant trajectories that were close to each other compared to an unweighted fitting procedure. Secondly, the introduction of the p-value provides us with a statistical quantity using which we may confirm or reject the null hypothesis.

4. The Orion nuclear cluster – a testbed

In order to carry out another test based on actual observations we applied our fitting procedure to other suggested runaways from the ONC which is close and young enough (McBride & Kounkel 2019) to study its runaway population. It has been studied quite extensively and Schoettler et al. (2020) used Gaia DR2 measurements to search for runaways and walkaway candidates from the ONC. Their list of stars were based on 2D velocity tracebacks where the ONC was kept stationary. In addition to the traceback, they also used isochrone measurements to shortlist candidates based on their age. A final addition was that of radial velocities which were taken from Gaia DR2 (Katz et al. 2019). Unlike the sample of massive early-type stars of Hoogerwerf et al. (2001), the sample consists of late-type stars. Consequently, there were 10 stars reported to be three-dimensional runaway candidates using these radial velocities. In addition, we supplemented the Gaia DR2 (Lindegren et al. 2018) with EDR3 measurements (Lindegren et al. 2021) along with Gaia DR2 radial velocities (as used by Schoettler et al. 2020). Similarly, Farias et al. (2020) provide a list of 16 stars of late spectral type also based on Gaia DR2 astrometry with corresponding Gaia DR2 radial velocities. Their method also relies on multiple flags like age and young stellar object color and magnitude, 2-d closest distance of approach, and three-dimensional overlap. Table 3 list p-values from our calculations for all these 26 stars based on both Gaia DR2 and EDR3 astrometry, where the lower limit was fixed to D_p = 5 pc that is the extent of the ONC. We also report the time and minimized distances for these stars. Of particular note here is that the stars have significantly high p-values. This result is consistent both with our method finding the minimized distance, and with the method used by Schoettler et al. (2020). In their paper, the mean values of the stars are used to find the velocities in the ONC rest frame. Since this method relies on the mean values to shortlist runaway stars it is consistent with high p-values since the observed parameters tend to be closer to their observed means.

The minimized times of flight we find are also within the same range as reported by the authors, and especially in the case of stars numbered 4, 6, and 24 the minimized times are much less than 1 Myr, meaning that they might indeed be within the cluster. However, for the latter star the minimized distance with EDR3 is 7.5 pc, which may hint at this being a nearby star and not a star from the cluster. Furthermore, while most of the stars have minimized distances below or near the nominal radius of 2.5 pc, the stars 12, 16, and 17 had close approaches greater than 40 pc, implying that these stars might not have their origin in the ONC. For star 15, the 20 pc difference between DR2 and EDR3 may be explained by the different parallaxes of 3.13 and 3.27 in DR2 and EDR3 respectively.

5. The Hoogerwerf et al. sample revisited

Hoogerwerf et al. (2001) in their seminal paper curated a list of 21 runaway stars based on HIPPARCOS astrometry along with their hypothesized places of origins. These origin sites were either nearby associations or open clusters, using the original HIPPARCOS catalog (ESA 1997) for parallaxes and proper motions of the stars. This catalog has been re-reduced by van Leeuwen (2007) and more precise measurements have been obtained by the Gaia satellite (Gaia Collaboration 2016), in particular the data release 2 (DR2 Gaia Collaboration 2018) and the early data release 3 (EDR3, Gaia Collaboration 2021). We compile the astrometric data from the new HIPPARCOS reduction as well as the two Gaia data releases in Table A.1. As quality indicators we use the goodness of fit parameter F2 for the HIPPARCOS data and the renormalized unit weight error (RUWE) for Gaia, requiring |F2| < 5 and RUWE < 1.4, respectively. Furthermore, since the Gaia brightness limit saturates at 3.6 mag, we do not use Gaia data for stars brighter than this. Therefore, for the stars ζ Oph, ζ Pup, and ι Ori, we stick to HIPPARCOS astrometry. The HIPPARCOS based distance has been reconfirmed by Howarth & van Leeuwen (2019). The RUWE parameters of Gaia EDR3 of two more stars, HIP 18614 and HIP 57669 exceed the limit and, therefore, their astrometry is taken from Gaia DR2. For the rest of the stars we stick to Gaia EDR3. Finally, we add the star HIP 82286 (HD 151397) which was studied by Blauuw to this list as it is a disk runaway. The spectro-photometric distance and radial velocity for this star have been provided by Dimpel (2021), and we use Gaia EDR3 for the astrometry.

Figure 3 compares Gaia DR2 and HIPPARCOS parallaxes for the runaway stars of Hoogerwerf et al. (2001), where corresponding Gaia parallaxes are available. For most of these stars, Gaia provides an order of magnitude better data, except for the brightest stars.

Fig. 3. Comparison of HIPPARCOS parallax values used by Hoogerwerf et al. (2001) with Gaia DR2 values.

Since Gaia’s catalog was released, hundreds of new clusters have been identified (Cantat-Gaudin & Anders 2020; Cantat-Gaudin et al. 2018; Gaia Collaboration 2018; Liu & Pang 2019). Apart from new clusters, the parameters of known clusters have been updated to more precise values (see Table A.2), as can also be seen in Fig. 4. Therefore, it is timely to revisit the Hoogerwerf et al. (2001) sample and check for the validity of the predicted birth clusters of the stars. The more so, because radial velocity studies have shown several stars of that sample to be radial velocity variable binaries (Jilinski et al. 2010).

Fig. 4. Comparison of parallax values used by Hoogerwerf et al. (2001) with Gaia DR2 based values for clusters mentioned in Hoogerwerf et al. (2001).

5.1. Radial velocities

Since radial velocities are one of the 6 parameters required for the kinematic analysis, they have to be updated as well. Gaia does not provide radial velocities for hot stars. Therefore, ground based spectroscopic studies are necessary. In order to identify close binaries showing radial velocity variations multi-epoch spectroscopy is required to pin down the orbit and in particular the systemic velocity. Jilinski et al. (2010) found five stars, previously suggested to originate from the Sco-Cen OB association, to be radial velocity variable and, hence, to be binaries. Therefore, we have to replace the radial velocities adopted by Hoogerwerf et al. (2001) by their system velocities (see Table 4). The space velocity of the stars (with respect to Galactic rotation) as calculated from the listed radial velocities, Gaia EDR3 astrometry, the circular velocity of the disk at solar distance and the LSR velocity of Model I of Irrgang et al. (2013) are also given in Table 4.

Concerning radial velocities of relevant stellar clusters data are unfortunately scarce (see Table A.2), for about a dozen cases they are based on one to four cluster members, only. In those cases the uncertainties may be larger than anticipated.

Furthermore, lower space velocities resulted from the new radial velocities and EDR3 measurements for many stars.

The binaries HIP 3881, HIP 38455, HIP 42038, HIP 46950, and HIP 76013 are likely walkaways too, rather than runaway stars, because their space velocities do not exceed 30 km s⁻¹. For the same reason three single stars (HIP 38518, HIP 48943, and HIP 49934) should be considered to be walkaways as well. The new space velocity of HIP 69491 even decreased to 9.30 km s⁻¹, meaning that the star may just be a very slow walkaway star. Hence, eleven stars of the sample are confirmed as runaways only, of which three have space velocities exceeding 100 km s⁻¹. A more appropriate criterion is the relative velocity with respect to the parent cluster, which can be calculated only when the cluster has been identified. They are listed in Table 5. Accordingly, the classification as walkaway or runaway star remains unchanged.

5.2. The cluster origins

We ran our p-value based distance minimization procedure for the runaway stars and the parent clusters suggested by Hoogerwerf et al. (2001). In Table B.1 we demonstrate that we can reproduce the results of Hoogerwerf et al. (2001), if we use the same input data. These results do not hold up when old radial velocities and new Gaia values are used for previously known single and binary stars as shown in Table B.2. Similarly, we show that making use of the Gaia data and more recent radial velocity measurements changes the conclusions in many cases considerably.

5.3. Identification of parent clusters and associations

Using our new method and equipped with accurate astrometry along with new radial velocities we searched for candidate clusters. We have used the best values available from the Gaia DR2, EDR3 (Lindegren et al. 2021), and the new reduction of HIPPARCOS (van Leeuwen 2007) and searched through more than 1000 clusters compiled from Kuhn et al. (2019), Gaia Collaboration (2018), Soubiran et al. (2018), Conrad et al. (2017), Cantat-Gaudin & Anders (2020), and Cantat-Gaudin et al. (2018).

Furthermore, important indicators of stellar ejection from the cluster and the ejection mechanism are the stellar age, the kinematic age, and the cluster age. We compiled age estimates of multiple clusters. The results are summarized in Table 5 and shall be discussed for each individual star in the following. To estimate errors on the time of flight from these clusters we use the approximation:

$$\begin{array}{c}\hfill T_{\text{rel}}=\frac{D_{\text{rel}}}{V_{\text{rel}}}\end{array}$$(5)

where D_rel and V_rel are the relative distance and the relative velocity between the star and the cluster. We sampled the radial velocities, proper motions, and parallaxes 10 000 times for ever system, converted the coordinates to a Galactocentric Cartesian system, and looked at the star and cluster system within the frame of reference of the cluster. This allowed us to derive an estimate of the standard deviation of the time of flight of the system, assuming that a straight line connects the two objects in this frame of reference.

Finally, while our simulations take into account Gaia correlations between the parallax, the proper motion component in right ascension, and that in declination from (Lindegren et al. 2018, 2021), the fitting procedure is not equipped for this yet, due to the fact that we do not sample the parameters but rather fit them. Using mock data we demonstrate in Appendix B that the impact of the correlations on the simulation results is negligible at the Gaia precision level. Nevertheless, for objects where correlations are more than ±0.3 between any two parameters we also provide simulations of trajectories to proof that the trajectories for these stars are not affected. Furthermore, none of our stars (for which a corresponding cluster was found) had a correlation factor between any two parameters beyond ±0.5.

6. The Hoogerwerf et al. sample revisited: Individual objects

We shall discuss the individual objects and their suggested parent clusters grouped into Walkaway stars (Sect. 6.1) and runaway stars (Sect. 6).

6.1. Walkaway stars and their suggested parent clusters

In the following we shall discuss the walkaway stars and their suggested parent parent clusters individually, arranged according to in increasing space velocity (see Table 4).

6.1.1. HIP 26241 (ι Ori)

HIP 26241 (ι Ori) is an excentric binary moving with a velocity of about 10 km s⁻¹ at almost right angles with respect to AE Aur/μ Col. The very different masses of the components hint at different ages of the components. N-body simulations by Gualandris et al. (2004) suggested that an exchange interaction occurred in an encounter between two low eccentricity binaries with comparable binding energy, which caused an older star to be swapped into the original ι Ori binary and the system to be ejected from the ONC. We are unable to confirm the place of origin (the ONC) as well as any encounter with the star AE Aurgiae. This difference occurs due to the change of parallax of ι Ori from 2.46 mas in the Hoogerwerf sample to 1.4 mas in the second reduction (van Leeuwen 2007)², AE Aur shifts from 1.74 mas to 2.57 mas in parallax, and the encounter no longer holds.

6.1.2. HIP 46950

HIP 46950 is a binary star, and therefore the original BSS scenario might not hold for it as pointed out by Jilinski et al. (2010). We find that it may have originated in the Gama Vel/Pozzo 1 cluster in the Vel OB2 association as a walkaway star, just like the walkaway star HIP 38518. This cluster is 13 Myr old, which is about the same time we find for the kinematic age of this star, making the dynamical ejection scenario more probable.

6.1.3. HIP 69491

The eclipsing binary HIP 69491 was suggested to come from the Upper Centaurus Lupus association (UCL) or Cep OB6. However, the high space velocity of 77 km s⁻¹ was puzzling for a binary system. Therefore, neither BSS nor DES could be confirmed. The findings of Bakış et al. (2008) for a radial velocity of v_rad = −10.3 km s⁻¹ resolved one of the problems with this star system. It is no longer a runaway candidate. As the third slowest star in our sample, the star barely qualifies as a walkaway. We find two major clusters as possible origin: NGC 6475 and C1439-661. The corresponding kinematic age is 16 Myr or 11 Myr, respectively³. Since NGC 6475 is almost 300 Myr old and C1439-661 is 65 Myr old (Bossini et al. 2019), only C1439-661 qualifies as a possible parent cluster. Accordingly, the most probable mechanism is BSS.

6.1.4. HIP 42038

The star HIP 42038 was suspected to come from either the Upper Centaurus Lupus association or IC 2391 as a BSS candidate. Due to it being a radial velocity variable, these origins had to be called in question. We show that the evidence for IC 2391 to be a cluster of origin is flimsy when the new radial velocity is used, as also found by Jilinski et al. (2010). Instead, we find three other clusters which may be the origin. These are Cl VDBH 99, NGC 2451A, and Cl Platais 9. The flight times are much higher than expected, being 34, 36, and 20 Myr respectively. The clusters themselves are older, with Cl VDBH 99 being 81 Myr, NGC 2451A being 44 Myr old, and Platais 9 being 78 Myr old Dib et al. (2018), Bossini et al. (2019). Due to the difference of 10 Myr between the age of the cluster and the flight time, and the high p-value of ≈0.96, NGC 2451A is the best candidate cluster for this star. However, in all three cases, the probability of HIP 42038 being a BSS is higher due to the fact that the kinematic age is much less than the age of the cluster. Such a scenario is possible if a third star in the system underwent a supernova. The resulting binary in this case would be a walkaway. This scenario was proposed by Hoogerwerf et al. (2001) to explain the evolution of HIP 69491, an eclipsing binary. We conclude that it may also apply for HIP 42038.

6.1.5. HIP 38455

Tetzlaff et al. (2011) claimed that the star HIP 38455 is not a runaway since it is a binary with a systemic radial velocity of 29.5 km s⁻¹ as opposed to −31 km s⁻¹ used by Hoogerwerf et al. (2001). We use the new radial velocity and confirm that it does not come from Collinder 135 as previously suggested. However, due to its low space velocity of ≈13 km s⁻¹, it might be a walkaway star which came close to the moving group Platais 3 about 18 Myr ago. The cluster is relatively old but its age is poorly known (208${}_{-41}^{+122}$ Myr, Bossini et al. 2019). Despite the age difference the cluster should not be discarded before a better age estimate becomes available.

6.1.6. HIP 48943

HIP 48943 is a B-type emission line star with a high projected rotational velocity of 190 km s⁻¹ (Shokry et al. 2018), which was thought to come from the Lower Centaurus Crux (LCC). Due to the new radial velocity, this does not hold true anymore. We find the cluster NGC 6405 intersected the trajectory of the star 25 Myr ago with a p-value of 0.37. The cluster is 48 Myr old (Liu & Pang 2019). The stellar age is not known, but the high projected rotational velocity of the star gives credence to the BSS. Its low relative velocity makes it a walkaway star.

6.1.7. HIP 3881 (ν And)

For the spectroscopic binary HIP 3881, we find the cluster Cl Alessi 20 as a possible site of origin. The estimated time of flight is 15 Myr based on EDR3 values for the star⁴. Richer et al. (2021) estimate the age of the cluster at 9 ± 4 Myr, while Kharchenko et al. (2013) give an estimate of 38 Myr. Therefore, it is possible that Cl Alessi 20 is the original cluster of the star. Since the star is a walkaway binary, it may have been created in a BSS. Depending on which age estimate for the cluster is selected, it could have a DES origin as well.

6.1.8. HIP 38518

While Gam Vel was the prime candidate for this star as was suggested by Hoogerwerf et al. (2001), we find the cluster Gulliver 9 to have also crossed paths with the star almost 8 Myr ago. This is a more significant approach than the one to Gam Vel when EDR3 is considered. However, Gulliver 9 is also part of the Vel OB2 region, meaning that the star was ejected from one of the regions of the Vel OB2 association. The original suggestion of BSS is the likely scenario in this case due to the fact that both the clusters are slightly older than the star. Its low relative velocity (see Table A.2) makes it a walkaway star.

6.1.9. HIP 76013

The star HIP 76013 is a walkaway Be type star with an He-SdO companion. According to the mass estimates for the He-SdO (Wang et al. 2021), the stars would be about 30 Myr old. We find no possible candidate clusters for this star using Gaia data. To check whether there may be a possible candidate with HIPPARCOS, we rerun the procedure and find that the star was close to the cluster M 67, 26 Myr ago with a p-value of 0.51. The probable origin of this binary star would be the DES if the stellar age were to be close to 26 Myr. However, the cluster is more than 3.5 billion years old (Ilin et al. 2021) and, therefore, it can not be a possible origin. The only other cluster we find with HIPPARCOS is the cluster Blanco 1, which was close 36 Myr ago with a p-value of 0.166. The cluster itself is about 100 Myr old, which may rule out an origin. However, due to the fact that the companion is an He-SdO, the visible star may have gone through a rejuvenation phase via mass transfer, and could be older than 30 Myr. A scenario of rejuvenation has been proposed before, for instance by Brown et al. (2010) to explain the origin of the hyper-velocity star HE 0437-5439.

6.1.10. HIP 49934

For the Be star HIP 49934, an origin in IC2602/IC2391 was suggested within the context of BSS. We find this origin does not hold, and find a different cluster NGC 3532, which intersected the star’s motion 5 Myr ago. However, the age of the cluster has been determined to be 400 Myr (Bossini et al. 2019), which is higher than our cut-off age. In addition, it has been found to host white dwarfs (Prišegen et al. 2021), cooling times seem to be consistent with the isochrone age. According to its low space velocity the star is a walkaway. The BSS does seem to be more probable in this case, as supported by the high projected rotational velocity of 220 km s⁻¹ (Balona & Ozuyar 2020) but the true cluster of origin is yet to be found.

6.1.11. HIP 57669

The Be star HIP 57669 was thought to originate in the Upper Scorpius region in the cluster IC 2602. Due to the age of the cluster, the flight time, and the high projected rotational velocity, the star was thought to be a BSS candidate, although no pulsar was traced back to it. The origin in IC 2602 was shown to not hold with new Gaia values (Neuhäuser et al. 2020a). We now find that its path crossed with the young moving group η Cha in the Chameleon star forming region 4–5 Myr ago. HIP 57669 travels slowly away from η Cha (20 km s⁻¹), which makes it a walkaway from that cluster. The age of the cluster is 11 ± 3 Myr (Dickson-Vandervelde et al. 2021; Bell et al. 2015), suggesting that the origin could indeed be a BSS origin. This would mean it was originally a triple system as Chini et al. (2012) found that HIP 57669 is a spectroscopic binary.

The star might also originate from the Cl Alessi 13 moving group, which is nearby the η Cha moving group, although the p-value is less than the 1σ limit. The flight time would be more than 10 Myr, compatible with the cluster age of about 40 Myr (Zuckerman et al. 2019). This would make HIP 57669 a BSS ejection, as originally suggested.

6.1.12. HIP 82868

For the single Be star HIP 82868, we find the moving group Alessi 13/X¹ Fornacius as a place of origin. This is a well studied cluster due to its proximity to Earth (≈100 pc), and new Gaia data has put the cluster age at 40 Myr (Zuckerman et al. 2019). The flight time of 12.5 Myr provides evidence for this star to be a BSS candidate. This is in agreement with the original suggestion, as well as the high projected rotational velocity of 370 km s⁻¹ (Jilinski et al. 2010). The star is a walkaway traveling at 27 km s⁻¹ from Alessi 13, and due to a high correlation the simulated trajectories are provided in Fig. 5. While this group may also be one of the candidates for HIP 57669, the two stars themselves never come close enough to each other with Gaia values.

Fig. 5. Histograms for the Monte Carlo simulations of back-traced trajectories for HIP 82868 and Alessi 13. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Alessi 13.

6.2. Runaway stars and their suggested parent clusters

In the following we shall discuss the runaway stars in the Hoogerwerf sample and their suggested parent cluster indiviudually arranged according to in increasing space velocity (see Table 4).

6.2.1. HIP 81377 (ζ Oph)

Hoogerwerf et al. (2001) suggested that the runaway star ζ Oph was ejected in a BSS event in the Upper Scorpius association and identified the pulsar PSR J1932+1059 as the runaway neutron star remnant of the former primary. A kick velocity of ≈350 km s⁻¹ to the neutron star disrupted the binary in the supernova explosion. Recently, however, Neuhäuser et al. (2020a) suggested that ζ Oph is linked to the pulsar B1706−16 in the association of stars known as the Scorpius-Centaurus-Lupus. This conclusion is supported by our calculations, although we find a slight mismatch of the best possible times, since a common origin with PSR J1932+1059 was found 1.3 Myr ago with a p-value of 0.41. However, the RV of the PSR is not well constrained. This makes it possible for the time of flight to be less, and therefore more cannot be said on the subject till then.

6.2.2. HIP 14514 (53 Ari)

The β Cep pulsator HIP 14514 is a classical runaway star (Blaauw 1961) and is confirmed to come from the Orion OB1 association, which it left almost 5 Myr ago. Due to the age of the cluster it is possible that 53 Ari was released as part of a binary system, but more data on the age of the star is required before this can be confirmed.

6.2.3. HIP 39429

We find that the O4 supergiant ζ Pup did indeed pass Trumpler 10, 2.5 Myr ago, which is consistent with its young age (3.2${}_{-0.2}^{+0.4}$ Myr Bouret et al. 2012).

Trumpler 10 is 55 Myr old (Bossini et al. 2019), which is more than the 30 Myr limit. However, Hoogerwerf et al. (2001) pointed out that ζ Pup may be a blue straggler, in which case a BSS origin in Trumpler 10 can not be ruled out.

6.2.4. HIP 102274

The B5 star HIP 102774 was thought to originate in the Cep OB2 association. We are unable to confirm this within 2σ certainty. We find two possible clusters, Alessi-Teutsch 5 and NGC 7160, from where the star might have originated. The two clusters are around 12.5 and 9.0 Myr old (Kounkel & Covey 2019; Dib et al. 2018; Conrad et al. 2017), meaning that the originally suggested BSS holds.

6.2.5. HIP 18614 (ξ Per)

HIP 18614 is a classical runaway star (Blaauw 1961), which is confirmed to come from the Per OB2 association via a DES ejection.

6.2.6. HIP 109556

We were unable to confirm the origin of the blue supergiant HIP 109556 (λ Cep) from the Cep OB3 association when using newer data from Conrad et al. (2017) and Dib et al. (2018). However, using data from Mel’Nik & Dambis (2009) allowed us to find a close encounter 2.95 Myr ago, which is less than that suggested by Hoogerwerf et al. (2001) and more in line with the expected age of the supergiant. However, Cep OB3b is less than 2.5 Myr old (Getman et al. 2018), meaning BSS would not hold. Therefore, Cep OB3a is a more likely parent.

6.2.7. HIP 22061 and HIP 29678

HIP 22061 and HIP 29678 are prime candidates for a dynamical ejection, from the cluster Collinder 69. However, the stars never come closer than 7 pc, which is incidentally also the radius within which 50% of the stars in Collinder 69 reside. Therefore, the cluster Collinder 69 remains the best candidate for these two stars. While Hoogerwerf et al. (2001) claimed the two stars were not related to the supernova that formed the ring around λ Ori, newer results have shown this explosion to be about 1 Myr ago (Dolan & Mathieu 2002). This is the same as the flight times of HIP 22061 and HIP 29678. Since the stars come from the slightly outer regions of the cluster, they may have been involved with this supernova.

For HIP 29678, Fig. 6 illustrates the histograms for close trajectories. Like our fitting procedure predicted, the simulations have almost no trajectories with closest distances lower than 5 pc.

Fig. 6. Histograms for the Monte Carlo simulations of back-traced trajectories for HIP 29678 and Collinder 69. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Collinder 69.

6.2.8. HIP 24575 (AE Aur), HIP 27204 (μ Col), and Collinder 69

The star μ Col (HIP 27204) has been one of the exemplary stars for the DES along with AE Aur (see Sect. 3.4). The two stars were found to move in almost opposite directions away from the ONC with velocities of about 100 km s⁻¹. However, we found evidence that μ Col might instead have its origin in Collinder 69, which is also in the Orion region, and thought to be where the two stars HIP 22061 and HIP 29678 originate from.

μ Col has a high p-value of 0.54 for a common encounter with the Collinder 69 cluster 3 Myr ago. The cluster is ≈5 Myr old (Kounkel 2020), and a prime site for the DES mechanism as evidenced by the two other stars. That μ Col can also be traced back to the ONC is a coincidence arising from the facts that the ONC is not that far away from Collinder 69 and AE Aurigae originated in the ONC at the same time. The fact that a runaway star can be a visitor to another dense star forming region has been shown before by Schoettler & Parker (2021), for the ONC itself. While they found past and future visitors to the ONC, μ Col was not considered, due to the fact that it has long been accepted in literature that it originated in the ONC⁵. Since the evidence for μ Col and its ONC origin have relied on Monte Carlo simulations, like the ones presented in Sect. 3.4, we also carried out MC simulations for Collinder 69 and μ Col to have been at the same place. The histogram of closest encounter from a million backward trajectories for each object, along with the time of closest encounter, and the ejection velocities, are shown in Fig. 7. The histograms match the spread for two objects thought to come from the same place (cf. Figs. 2 and B.1), with the peak of the distribution below 20 pc. Since this common encounter is older than 2.5 Myr, it is more probable for μ Col to have been born in Collinder 69 instead of the ONC. The best-fit trajectories for both the stars and their respective places are illustrated in Fig. 8. The average time of more than 3 Myr provides evidence that μ Col may have originated in this cluster instead of ONC.

Fig. 7. Histograms for the Monte Carlo simulations of back-traced trajectories for μ Col and Collinder 69. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Collinder 69.

Fig. 8. Orbits for AE Aurigae and μ Col using best-fit parameters for their respective origins from the ONC and Collinder 69. The figure illustrates how μ Col may have been ejected from Collinder 69, and while it passed through the ONC, AE Aurigae was ejected from the latter. X and Y represent the Galactocentric Cartesian coordinates with the Sun at (8.4 kpc, 0). Crosses represent 0.35 Myr of backward trajectory. The trajectories begin at the respective labels. Trajectories are plotted using Galpy (Bovy 2015).

6.2.9. HIP 91599

The star is a B-type star (T_eff = 26 530 K, log g = 3.65 Daflon et al. 2004) with a projected rotational velocity of 92 km s⁻¹ (Daflon et al. 2007; Gonneau et al. 2020). Its parallax given in the HIPPARCOS catalog and its new reduction (ESA 1997) is large (3.61 ± 1.16 mas and 3.41 ± 0.72 mas, respectively). However, its parallax from both Gaia catalogs are much smaller (0.93 mas and 0.91 mas in Gaia DR2 and EDR3, respectively). The star is in a crowded field at a very low Galactic latitude of −1.5°. This allows an independent distance estimate be made from reddening-distance relations. By fitting the spectral energy distribution (SED) and color indices of the star (see Fig. 9), we derived its interstellar reddening as described in Heber et al. (2018). Our fit gives a color excess of E(B − V) = 0.77 ± 0.02 mag. The reddening-distance relations of Capitanio et al. (2017) suggest that its distance must exceed 1 kpc consistent with the Gaia parallax. Consequently, its space space velocity is high, making HIP 91599 the second fastest star in the disk runaway sample.

Fig. 9. Comparison of synthetic and observed photometry for HIP 91599. The top panel displays the spectral energy distribution. The colored data points are filter-averaged fluxes that were converted from observed magnitudes. The dashed horizontal lines depict the full widths at 10% of the maximum of the filter transmissions, while the solid gray line represents the best fitting model (spectral resolution of 6 Å). In order to reduce the steep slope of the SED the flux is multiplied by the wavelength to the power of three. Bottom panel: residuals, χ: differences between synthetic and observed magnitudes divided by the corresponding uncertainties. The photometric systems have the following color code: UV magnitudes from the TD1 satellite (violet; Nandy et al. 1978) labeled with central wavelengths, optical magnitudes from Tycho2 (Høg et al. 2000), HIPPARCOS (van Leeuwen 2007), Gaia (blue Gaia Collaboration 2020), and Mermilliod (2006); IR photometry from 2MASS (red Cutri et al. 2003) and DENIS (yellow Denis 2005). Right hand panel: color indices in the Johnson (cyan) and Geneva (red) systems (Mermilliod 2006; Rufener 1988).

HIP 91599 was thought to come from either the Per OB2 or the Per OB3 associations. We find that the cluster NGC 6883 has a minimized orbit with the star, 10.7 Myr ago with a p-value of 0.92 (we adopt a distance of 1380 pc from Dib et al. (2018), assuming a 10% error). This cluster is around 50 Myr old (Dib et al. 2018; Conrad et al. 2017), which is more than a 30 Myr difference. However, HIP 91599 is another blue straggler candidate (Hoogerwerf et al. 2001), therefore making BSS a likely origin mechanism from this cluster despite the time difference. This is further supported by the high Helium abundance of the star. The best-fit orbit is shown in Fig. 10.

Fig. 10. The orbits of the star HIP 91599 and the cluster NGC 6883 for the best-fit parameters shown here for 10.7 Myr in the past. Crosses represent 1.1 Myr of backward trajectory and circles represent present day positions. X and Y are Galacto-centric Cartesian coordinates. Trajectories are plotted using Galpy (Bovy 2015).

6.2.10. HIP 82286 (HD 151397) – an overlooked runaway B star

The star HIP 82286 (HD 151397) was part of the runaways already studied by Blaauw (1961), who suggested the cluster NGC 6231 as a possible site of origin. The star was not included in the (Hoogerwerf et al. 2001) sample, because HIPPARCOS failed to deliver useful astrometry. The new reduction improved it, but the uncertainty of the parallax still remained at the 70% level. Its high radial velocity of v_rad = 151 ± 4.7 km s⁻¹ (de Bruijne & Eilers 2012) makes it a highly interesting object. Hence the star was revisited by Dimpel (2021) in a spectroscopic analysis of high Galactic latitude runaway stars (see Sect. 7), who found it to be a very young (1.3 Myr), hot (T_eff = 29 760 ± 300 K) and massive (13.8 M_⊙) main sequence star.

NGC 6231 was not confirmed as its parent cluster, and since then no possible cluster or association has been linked to the origin of this star. We used the radial velocity along with Gaia astrometry to search for possible origin sites for this star. About 0.7 Myr ago, HD 151397 was within 7 pc of the cluster Cl VDBH 205 (in the SCO OB1 association) with a high p-value of more than 0.96 when using Gaia EDR3 astrometry. This is slightly outside the inner part of the cluster which has 50% stars within 3 pc. The kinematic age is less than the star’s evolutionary age of 1.5 Myr. The cluster is somewhat older. Its age is given as 13 Myr by Sampedro et al. (2017) and 6.8 Myr (Tarricq et al. 2021), both, however, without uncertainties. It is possible that the star is a blue straggler resulting from mass transfer in a binary system prior to the SN explosion. Since this star has a correlation higher than 0.3 between the proper motions, we show the simulations for close encounters in Fig. 11. Therefore, we suggest that this star may have originated in Cl VDBH 205 in the BSS.

Fig. 11. Simulated histograms of closest encounter, time of closest encounter, and ejection velocities for HD151397 and Cl VDBH 205 for the fraction of trajectories coming closer than 20 pc.

6.3. The Hoogerwerf et al. (2001) sample in the Gaia era

The results described above for the stars of the Hoogerwerf et al. (2001) individually are summarized here. Firstly, the space velocities of several stars were found to be lower than previously derived, leaving us with a sample of 12 walkaway with space velocities between 5 and 30 km s⁻¹ and 11 runaway stars with space velocities between 37 and 120 km s⁻¹. Furthermore, there are eight spectroscopic binaries among the walkaways and only one (ξ Per) among the runaways. Parent clusters or associations have been identified for all but 4 runaway stars (HIP 76013, HIP 109556, HIP 38455, and HIP 49934). For HIP 38455 and HIP 49934 we found spatial coincidence in the past with two clusters, but the cluster ages ruled out a common origin. For HIP 109556, the only close encounter with Cep OB3 was found when using HIPPARCOS second reduction based data from Mel’Nik & Dambis (2009). We confirm the originally suggested parent clusters for 53 Ari, HIP 22061 and HIP 29678, AE Aurigae, ψ Per, and HIP 38518.

We found that the ι Ori was not involved in a close encounter with AE Aurigae or the ONC. Furthermore, μ Col was found to have a close encounter with Collinder 69, 3 Myr ago, earlier than its flight time from the ONC. Therefore, the classical example of dynamical ejection scenario is no longer found to hold true.

New parent clusters were identified for 11 stars. Out of these, 9 clusters had p-values higher than the 1σ limit. The other 2 had p-values between 1.2−1.5σ levels. Finally, ζ Oph was confirmed to be associated with the runaway neutron star PSR B1706-16 as has recently been suggested by Neuhäuser et al. (2020a).

7. Runaway B stars in the Galactic halo

The list of stars we have so far discussed are all stars with space velocities between 10 and 100 km s⁻¹, with a few slightly above 100 km s⁻¹. These stars are young and possibly still within the disk of the Galaxy. Many are Be stars, binaries, and supergiants. However, runaway stars far from the Galactic plane are also known (Greenstein et al. 1956; Greenstein & Sargent 1974; Tobin 1987). Silva & Napiwotzki (2011) studied the sample known at the time. In addition several runaway B-type stars have been discovered and analyzed since then (e.g., McEvoy et al. 2017; Irrgang et al. 2019; Li et al. 2018). Recently, a subset of main sequence B stars from that sample has been analyzed by Dimpel (2021) using high resolution spectra and Gaia based kinematics. All stars are B main sequence stars and none is known to be a binary. We refer to his results as listed in Table A.3 or to Dimpel (2021). Due to their high velocity, these stars are farther away than the runaways discussed so far, making it difficult to trace back their orbits. One of the most important use of our fitting procedure is to find possible sites of origins of the halo runaway stars, which is presented in this section. The kinematic parameters of the stars are summarized in Table A.3. Due to the fact that these stars are much farther away from us (> 1 Kpc), we use Gaia EDR3 proper motions, but stick to spectroscopic distances from Dimpel (2021).

7.1. PG 1332+137 and Sai 132

The B5V star PG 1332+137, also known as Feige 84, is a 5.5 M_⊙ star with an evolutionary age of 46 Myr and a Galactic rest frame velocity of 299 km s⁻¹ (Dimpel 2021). We find only one cluster close to it. The star rotates moderately fast (140 ± 15 km s⁻¹Lynn et al. 2004). The cluster Sai 132 was close to the star with a p-value of 0.99, 20.7 Myr ago. However, the cluster age is poorly known. According to Kharchenko et al. (2013) the cluster is old (≈500 Myr), while Liu & Pang (2019) give an age of just 6 ± 0.4 Myr. Hence, any conclusion has to await the resolution of the cluster age discrepancy.

7.2. PHL 159 and Cl Basel 2

The star PHL 159 is an 8 M_⊙ star with an evolutionary age of 21 Myr. We find the cluster Basel 2 was nearby the star with high p-value of 0.73. The time of flight of 24.9 Myr is slightly inconsistent with the evolutionary age, a discrepancy which has been pointed out before by Dimpel (2021). This tension may be accounted for if PHL 159 underwent a BSS ejection, which would make the star appear younger due to binary evolution. The cluster Basel 2 however is much older than the star (Glushkova et al. 2010), making it inappropriate for such a mechanism.

7.3. PB 5418 and Cl Gulliver 43

PB 5418 is a 6.5 M_⊙ star with an evolutionary age of 42 Myr old. It has a radial velocity of v_rad = 146.6 ± 2.6 km s⁻¹. We find that it crossed paths with the open cluster Gulliver 43 20 Myr ago with a p-value of 0.67. The age of the cluster is not known but the time of flight and evolutionary age suggest that BSS is the more probable scenario in this case.

7.4. HIP 70275 and NGC 1502

The 9 M_⊙ star HIP 70275 is ≈17 Myr old and has a radial velocity of v_rad = 242.3 ± 0.2 km s⁻¹. The star comes close to the young cluster NGC 1502 with a p-value of 0.52, 11.2 Myr ago. The age of the cluster is also around 11 ± 0.2 Myr (Yu et al. 2018). Along with the low projected rotational velocity (v_rot = 64.3 ± 3 km s⁻¹McEvoy et al. 2017) a DES ejection would be most probable if the cluster is indeed the correct one. Otherwise the difference in evolutionary age and flight time suggests a BSS ejection, which would not have been possible in NGC 1502, again assuming the cluster age is correct⁶.

7.5. HIP 11809 and [FSR2007] 0883

HIP 11809 (Feige 23) was close to the cluster [FSR2007] 0883, 10.2 Myr ago respectively. The estimated stellar age of 100 Myr, and the high projected rotational velocity of 240 km s⁻¹ (Martin 2004) implies a BSS origin. The cluster is older and the published ages range from 234 Myr (Dias et al. 2014) to 700 Myr (Liu & Pang 2019). Therefore it may not be the origin cluster of the star, if the cluster age estimates are correct.

7.6. HIP 114569 and Cl VDBH 217

The 5.5 M_⊙ star HIP 114569 was 25 pc away from the cluster VDBH 217, 7 Myr ago. This is much more than the radius of the cluster, but no other cluster was found to be as close to the star. Furthermore, the cluster is slightly younger than the star at 40 ± 10 Myr (Clariá et al. 2019; Sampedro et al. 2017), and has only one star with known Gaia RV. Compared to the evolutionary age of the star and its flight time, BSS is the most probable scenario, which further makes this cluster unlikely as the parent. Therefore, it is likely that we have not yet identified a cluster for this star.

7.7. HIP 105912 and Teutsch 145

The 10 M_⊙ star HIP 105912 had a close encounter 14 Myr ago with the older cluster Teutsch 145 with a p-value of 0.66. The star is 20 Myr old, leaving some time for binary evolution, and therefore BSS may be the origin here. However, firstly the radial velocity of the cluster relies on only one star, and secondly, age estimates for the cluster range from 200 Myr (Sampedro et al. 2017) to more than 500 Myr (Kharchenko et al. 2013; Joshi et al. 2016). Therefore, Teutsch 145 is unlikely to be the possible place of origin, if the cluster’s radial velocity and age estimates are correct.

7.8. BD-02 3766 and Kronberger 38

The 10 M_⊙ star was near the cluster Kronberger 38 with a p-value of 0.54, almost 15 Myr ago. This is also the estimated evolutionary age of the star, therefore making it a prime candidate for DES. Not much is known about the cluster itself, since it is in a region obscured by dust, and also only has one measured RV star.

7.9. BD-15 115 and Juchert 3/NGC 6755

The 8 M_⊙ star B-15 115 is 32 Myr old, and around the same time it was near 3 clusters: Juchert 3, NGC 6705, and NGC 6755 with p-values of 0.52, 0.33, and 0.01. However, NGC 6705 is older than 300 Myr (Rain et al. 2021; Liu & Pang 2019) and therefore not a suitable choice for DES. On the other hand, not much is known about the cluster Juchert 3, except that is a slightly younger cluster (Kharchenko et al. 2013). However, Juchert 3 has only 1 RV star measured (Soubiran et al. 2018). Finally, the lowest p-value cluster NGC 6755 is around 174 Myr old (Liu & Pang 2019). Therefore, only Juchert 3 is a good candidate.

8. The hyper-runaway candidate HIP 60350 and Cl Pismis 20

The star HIP 60350 was considered as a hyper-runaway candidate with a space velocity of about 400 km s⁻¹ by Irrgang et al. (2010). It was discussed by Tenjes et al. (2001) as having a possible origin in the cluster NGC 3603. However, due to the difference in the cluster age and the stellar age, the evidence for this origin was quite weak and was disproven by Irrgang et al. (2010). In their work, Irrgang et al. (2010) presented a detailed kinematic and spectroscopic analysis of the star and suggested clusters which were close to the location of the star when it had crossed the disk of the Milky Way. However, a definitive place of origin is still lacking, although EDR3 data has constrained the stellar kinematics (Irrgang et al. 2021). We adopt a radial velocity of v_rad = 263.6 ± 0.5 km s⁻¹ as used by Irrgang et al. (2021), and use the proper motions according to Cantat-Gaudin & Brandt (2021).

We used our fitting procedure and were able to falsify all previously suggested origins in those clusters older than the stellar age (Irrgang et al. 2021 place the star at a stellar age of 50 Myr, much older than the previously suggested clusters). The only cluster left is the cluster NGC 5617, for which we get a p-value of 0.002 for a close encounter of less than 1 pc, which while sufficient might not be the best. Similar to Sect. 4, we searched again for a new cluster using Gaia EDR3 values, and found close encounters of the star with the clusters Cl Pismis 19 and Cl Pismis 20, with p-values of 0.40 and 0.43 (see Table 6). The first encounter happened 14.3 Myr ago while the second was 14.98 Myr ago, same as the time of disk crossing of 14.9 ± 1.0 suggested by Irrgang et al. (2010). Cl Pismis 20 is slightly younger than 50 Myr (Dib et al. 2018), but may still be a suitable choice for the origin⁷ On the other hand Cl Pismis 19 is a much older cluster of around 1 Gyr of age (Angelo et al. 2021) and is therefore not a viable candidate.

Since the stellar age is much older than the kinematic age, this star was possibly involved in a supernova explosion. According to Irrgang et al. (2010), one explanation could be the involvement of a high mass star which would be a Wolf-Rayet star of stellar type WN. Such stars have been known to be hosted by clusters, such as Cl Pismis 20 (Rate et al. 2020). The mass required for the Wolf-Rayet star would be greater than 15 M_⊙. Such a star would now be either a neutron star or a black hole, rendering it difficult to find the original companion of HIP 60350. Figure 12 shows the orbit of the star and the clusters, illustrating the high latitude nature of the star and how it may have been ejected from the disk. The results for HIP 60350 along with all the other halo runaways are summarized in Table 6.

Fig. 12. The 3 dimensional orbit of the star HIP 60350 shown here in X–Y and X–Z Cartesian coordinates. The two clusters Cl Pismis 19 and Cl Pismis 20 are also plotted for 14.3 and 14.9 Myr respectively. Crosses represent 1.5 Myr in time while circles represent the present day positions. Trajectories are plotted using the revised halo potential of Model I from Irrgang et al. (2013) in Galpy (Bovy 2015).

9. Summary and conclusion

The excellent quality of Gaia astrometry opens up the path to study the kinematics of young runaway stars to such a high precision that the place of origin in open stellar clusters can be identified uniquely even when it is a few kiloparsec away. Moreover, hundreds of new clusters have been discovered utilizing Gaia data.

To cope with the quality and growth of observational data, we developed an efficient minimization method to calculate whether two or more objects may come from the same place. Our method is based on a fitting procedure which fits the trajectories of objects by minimization of the distance between them taking into account the significance of the trajectories by means of the statistical p-value. The 2σ limit for 8 parameters is set, and any p-value below this level allows us to reject the null hypothesis and conclude that the two or more objects do not come from the same place with 2σ certainty. We test our procedure first on the classical runaways AE Aurigae and μ Columbae thought to have originated in the ONC 2.5 Myr ago. Furthermore we confirm the results of Schoettler et al. (2020) and Farias et al. (2020) for stars thought to originate in the ONC, except for three stars that we found to have closest approaches of more than 40 pc from the ONC center.

So far, the most sophisticated study of cluster runaways was carried out by Hoogerwerf et al. (2001) based on HIPPARCOS data. Using new data from Gaia DR2 and EDR3, and HIPPARCOS new reduction, along with new radial velocities, we run our procedure on the runaway sample of Hoogerwerf et al. (2001). We confirm the previous results for 53 Ari, AE Aurigae, ψ Per, HIP 22061, and HIP 38518. Among the rest of the sample we find new candidate clusters as sites of origin using new Gaia data for most of them. In many cases the reason for the new identifications can be traced back to the radial velocities used in the original study. Follow-up spectroscopy had shown that some stars are radial velocity variable binaries. We used the systemic velocities, which had a strong impact on the kinematic results. Only half of the sample is confirmed to be runaway stars, while the space velocities of the second half are so low (< 30 km s⁻¹), that they do not qualify as runaways, but have to be considered walkaway stars. In particular we show that the classical runaways AE Aurigae and μ Columbae might not have originated together, with μ Columbae having an earlier ejection from Collinder 69, a cluster near the ONC. Furthermore, the binary system ι Ori is found to not have had a close encounter with either AE Aurigae or the ONC. We use cluster age and flight time of the stars to investigate whether the ejection was due to a binary supernova or due to a dynamical ejection.

We use our procedure on a list of runaway B stars in the halo compiled from Silva & Napiwotzki (2011) using results from recent analyzes by Dimpel (2021). We are able to find candidate clusters for four stars where these stars and derive their ejection velocities to lie in the range 160–410 km s⁻¹ (Dimpel 2021) considerably larger than for the Hoogerwerf et al. (2001) sample. Three other star-cluster pairs may have a common origin, but due to nonavailability of ages for these clusters, we are unable to confirm this. We also find that in some cases the clusters may have been at the same places as the stars in the past, but the cluster ages rule out any common origin.

Finally, the hyper-runaway B-star HIP 60350 was revisited. The recent Gaia EDR3 based kinematic analysis by Irrgang et al. (2021) already dismissed the star as hyper-runaway star because its Galactic restframe velocity is lower than the local escape velocity. All of the previously proposed parent clusters had to be dismissed, and based on both spatial and temporal coincidence, Pismis 20 turn out to be the only good candidate. The ejection velocity of 338 km s⁻¹ is lower than previously derived but is similar to that of the Galactic Halo runaway B stars.

While radial velocities of the runaway stars are very well known, this is not the case for many open clusters, whose radial velocities are based on very few stars, occasionally on one star, only. This is especially true for open clusters more than 1 kpc away from the Sun. Additional radial velocity surveys of open clusters are needed to derive accurate cluster radial velocities, which would further allow their kinematics to be constrained better.

Acknowledgments

A.B., A.I. and U.H. acknowledge funding by the Deutsche Forschungsgemeinschaft under grants HE 1356/70-1 and IR 190/1-1. We would like to thank the anonymous referee for their insightful comments. 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 has made use of the SIMBAD database, operated at CDS, the VizieR catalog access tool, CDS, Strasbourg, France (DOI : 10.26093/cds/vizier), and TOPCAT table processing software (Taylor 2005).

References

Appendix A: Kinematic parameters for runaway stars, and open clusters used in this study

Table A.1 list the astrometric data from Hipparcos and Gaia for the sample of Hoogerwerf et al. (2001) (see Sect. 5), while Table A.2 compiles the data for the relevant open cluster. In Table A.3 the astrometric data for the runaway B stars are listed (see Sect. 7 and Sect. 8).

Appendix B: Testing the p-value method for the Hoogerwerf+ sample and mock data

For the purpose of testing our methods as well as to see how any predictions deal with uncertainties we created some mock data. For this we recreated an interaction between AE Aurigae and μ Columbae in the present time. We took the present position of the Orion Nebula cluster (ONC, Trapezium cluster) and assumed the two stars to be in the center of the ONC. We provided the same parallax to both the stars assuming they come from the same point in space. Then we applied our orbit calculator on the two stars using their present day velocities (this may seem nonphysical at first but the initial velocities are not as constrained in the interaction as the three-dimensional position). We ran the two orbits for 2.5 Myr (approximately the flight time of present day AE Aurigae and Mu Col.) and noted the positions and velocities of the two stars now. The common starting spatial parameters used for the two stars were:

α: 5h 35m 16s; δ : -5° 23’ 40”; ϖ: 2.50 mas *

Therefore, we ended up with two different trajectories and used the end points of these trajectories as the values for our mock data. To simulate Hipparcos and Gaia precision we add to all the final velocities and parallaxe in the first simulation 10% errors, which is typical for Hipparcos precision. In a second simulation we the same relative errors as in EDR3 catalog. Finally, we trace back the two stars in time again for both adopted precisions. To check whether our mock data are consistent, we apply the simulation technique by computing 100000 total backward trajectories and plot histograms shown in Fig. B.1. The results of the simulation using Gaia precision (lower panel of Fig. B.1) are close to that derived from the new observations shown in Fig. 2. The exercise demonstrates the great improvement to be achieved from high precision measurements.

Fig. B.1. Histograms for simulated back-traced mock orbits of two stars, with 10% errors (top) and EDR3 based relative errors on AE Aurigae and μ Col (bottom). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the stars with respect to each other.

Our fitting procedure found a p-value of 0.85, even in the simulation with 10% errors, close to the values found for real data of Hipparcos.

Our MC simulations allow us to include astrometric correlations, which are available for Gaia DR2 (Lindegren et al. 2018) and EDR3 (Lindegren et al. 2021) data. To test for the effects of such correlations, we added correlation factors of 0.9 between the parallax, proper motion in right ascension, and proper motion in declination. The histograms in Fig. B.2 show the effect of these correlations. For EDR3 based errors, the errors are small enough that the correlations do not lead to significant changes in the inference of the trajectories. A similar result is found when the correlation is 0.5.

Fig. B.2. Histograms for simulated back-traced mock orbits of two stars with EDR3 based relative errors on AE Aurigae and μ Col (bottom), when a correlation factor of 0.9 is added to all the three major parameters (top). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the stars with respect to each other.

Finally, we also ran our p-value based distance minimization procedure for the runaway stars and clusters to identify their places of origin suggested by Hoogerwerf et al. (2001) using the original radial velocities and proper motions in most cases as a proof of test for our method to compare our results with the published one. For clusters we utilize values from the Simbad database, based on new Gaia DR2 kinematics where available. Table B.1 shows our results.

All Tables

All Figures

	Fig. 1. Cluster extent.
In the text

Fig. 2. Histograms for back-traced orbits of AE Aurigae and μ Col from Monte Carlo simulations. with HIPPARCOS second reduction (top) and Gaia EDR3 data (bottom). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative velocity of the stars with respect to each other.

In the text

	Fig. 3. Comparison of HIPPARCOS parallax values used by Hoogerwerf et al. (2001) with Gaia DR2 values.
In the text

	Fig. 4. Comparison of parallax values used by Hoogerwerf et al. (2001) with Gaia DR2 based values for clusters mentioned in Hoogerwerf et al. (2001).
In the text

Fig. 5. Histograms for the Monte Carlo simulations of back-traced trajectories for HIP 82868 and Alessi 13. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Alessi 13.

In the text

Fig. 6. Histograms for the Monte Carlo simulations of back-traced trajectories for HIP 29678 and Collinder 69. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Collinder 69.

In the text

Fig. 7. Histograms for the Monte Carlo simulations of back-traced trajectories for μ Col and Collinder 69. The left most panel shows the fraction of distances of closest encounters between the star and the cluster. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the star with respect to Collinder 69.

In the text

Fig. 8. Orbits for AE Aurigae and μ Col using best-fit parameters for their respective origins from the ONC and Collinder 69. The figure illustrates how μ Col may have been ejected from Collinder 69, and while it passed through the ONC, AE Aurigae was ejected from the latter. X and Y represent the Galactocentric Cartesian coordinates with the Sun at (8.4 kpc, 0). Crosses represent 0.35 Myr of backward trajectory. The trajectories begin at the respective labels. Trajectories are plotted using Galpy (Bovy 2015).

In the text

Fig. 9. Comparison of synthetic and observed photometry for HIP 91599. The top panel displays the spectral energy distribution. The colored data points are filter-averaged fluxes that were converted from observed magnitudes. The dashed horizontal lines depict the full widths at 10% of the maximum of the filter transmissions, while the solid gray line represents the best fitting model (spectral resolution of 6 Å). In order to reduce the steep slope of the SED the flux is multiplied by the wavelength to the power of three. Bottom panel: residuals, χ: differences between synthetic and observed magnitudes divided by the corresponding uncertainties. The photometric systems have the following color code: UV magnitudes from the TD1 satellite (violet; Nandy et al. 1978) labeled with central wavelengths, optical magnitudes from Tycho2 (Høg et al. 2000), HIPPARCOS (van Leeuwen 2007), Gaia (blue Gaia Collaboration 2020), and Mermilliod (2006); IR photometry from 2MASS (red Cutri et al. 2003) and DENIS (yellow Denis 2005). Right hand panel: color indices in the Johnson (cyan) and Geneva (red) systems (Mermilliod 2006; Rufener 1988).

In the text

	Fig. 10. The orbits of the star HIP 91599 and the cluster NGC 6883 for the best-fit parameters shown here for 10.7 Myr in the past. Crosses represent 1.1 Myr of backward trajectory and circles represent present day positions. X and Y are Galacto-centric Cartesian coordinates. Trajectories are plotted using Galpy (Bovy 2015).
In the text

	Fig. 11. Simulated histograms of closest encounter, time of closest encounter, and ejection velocities for HD151397 and Cl VDBH 205 for the fraction of trajectories coming closer than 20 pc.
In the text

Fig. 12. The 3 dimensional orbit of the star HIP 60350 shown here in X–Y and X–Z Cartesian coordinates. The two clusters Cl Pismis 19 and Cl Pismis 20 are also plotted for 14.3 and 14.9 Myr respectively. Crosses represent 1.5 Myr in time while circles represent the present day positions. Trajectories are plotted using the revised halo potential of Model I from Irrgang et al. (2013) in Galpy (Bovy 2015).

In the text

Fig. B.1. Histograms for simulated back-traced mock orbits of two stars, with 10% errors (top) and EDR3 based relative errors on AE Aurigae and μ Col (bottom). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the stars with respect to each other.

In the text

Fig. B.2. Histograms for simulated back-traced mock orbits of two stars with EDR3 based relative errors on AE Aurigae and μ Col (bottom), when a correlation factor of 0.9 is added to all the three major parameters (top). The left most panel shows the fraction of distances of closest encounters between the two stars. The middle panel is the histogram of the times for the closest encounters of the two stars, and the right most panel plots the relative ejection velocity of the stars with respect to each other.

In the text