© The Authors 2023

1. Introduction

Massive early-type OB stars are the most luminous stars in the Milky Way. They are relevant for the study of the kinematics and dynamics of the young stellar populations, the metallicity enrichment in the Galaxy, and the feedback in the interstellar medium through gamma-ray bursts and supernova explosions, among many other topics (see, e.g., Vanbeveren et al. 1998; Woosley & Bloom 2006). O stars are characterized by powerful stellar winds and short life-spans (< 10⁷ yr), while B stars have longer life-spans (up to ∼10⁸ yr). Be stars are a particular type of B stars that display Balmer emission lines and an excess of infrared emission because they present circumstellar envelopes in the form of decretion disks (Slettebak 1988; Rivinius et al. 2013). A relevant feature of OB stars is that they are mostly found in binaries, while at least 70% of them form interacting binaries (Chini et al. 2012; Sana et al. 2012; Moe & Di Stefano 2017). After stellar evolution, these systems can form high-mass X-ray binaries (Lewin & van der Klis 2006), gamma-ray binaries (Dubus 2013), millisecond pulsar systems, and double neutron stars (van den Heuvel 2007), which allow the study of nonthermal processes.

On the other hand, more than 30% of the O stars and about 5–10% of the B stars are considered to be runaway stars. These stars have a high peculiar velocity with respect to their environment (Blaauw 1961; Stone 1979; Boubert & Evans 2018). There are two possible scenarios to explain the existence of runaway stars: (1) the dynamical ejection scenario (DES), where a close three- or four-body interaction in the core of a dense cluster ejects a massive star (Poveda et al. 1967), and (2) the binary supernova scenario (BSS), where the mass loss in the supernova (SN) explosion of the most evolved star in a binary system induces either a runaway velocity to the remaining binary system, which now contains a compact object, or disrupts the binary system and induces runaway velocities to both the remaining massive OB star and the remaining compact object (Blaauw 1961). Later works explored the idea that massive runaway stars could host compact companions (Bekenstein & Bowers 1974; Stone 1982; van Oijen 1989). There is now observational evidence of runaway X-ray binaries (see, e.g., van den Heuvel et al. 2000). In addition, three out of the ten known gamma-ray binaries are runaways: LS 5039, PSR B1259−63, and 1FGL J1018.6−5856 (Marcote et al. 2018 and references therein). Therefore, runaway massive stars can be part of close binary systems. In a sample of O runaway stars, Chini et al. (2012) obtained a binary fraction of 69%.

Several works that searched for massive runaway stars can be found in the literature, from early efforts based on radial and space velocities from the literature (Cruz-González et al. 1974; Stone 1979) and those using HIPPARCOS proper motions (Hoogerwerf et al. 2001; Mdzinarishvili 2004; de Wit et al. 2005; Tetzlaff et al. 2011) to recent publications using Gaia DR1, DR2, and EDR3 data (Maíz Apellániz et al. 2018; Boubert & Evans 2018; Kobulnicky & Chick 2022). Not only is there diversity in the input data, but also in the runaway detection method. Some authors based the criterion for detecting runaways on the three-dimensional (3D) peculiar velocity of the stars above a velocity threshold, which is usually in a range from 25 to 40 km s⁻¹ (Blaauw 1961; Cruz-González et al. 1974; Hoogerwerf et al. 2001). Other authors complemented this criterion with that of the distance of the stars above or below the Galactic plane, setting |Z|> 0.25–0.50 kpc as the limit (de Wit et al. 2005; Tetzlaff et al. 2011). Because precise radial velocities and distances are lacking, other authors used a two-dimensional (2D) proper-motion-only method, which allows detecting about half of the existing runaways (Maíz Apellániz et al. 2018). Kobulnicky & Chick (2022) used a method in which a star is considered as runaway when it has a 2D peculiar velocity above 25 km s⁻¹ with respect to its environment, assuming a Galactic rotation curve. This diversity of input data and methods translates into a diversity of results.

Our aim in this work is to search for new massive O and Be runaway stars, some of which could be (new) high-mass X-ray or gamma-ray binaries. In particular, the discovery of new gamma-ray binaries may shed light on some of the unresolved questions in these systems (Dubus 2013; Dubus et al. 2017). To this end, we started to search for runaway massive stars using O-star catalogs and Gaia DR2 data using a 2D velocity method (Ayan-Míguez & Ribó 2019), as well as O- and Be-star catalogs and Gaia EDR3 data (Carretero-Castrillo et al. 2023). In this work, we continue this project, but use Gaia DR3 data and improve several aspects of the method.

This paper is organized as follows. In Sect. 2 we describe the different catalogs we used. In Sect. 3 we explain the cross match of the massive star catalogs with Gaia DR3, the quality cuts we applied to the data, and the properties of the resulting catalogs. In Sect. 4 we describe our method, how we computed the velocities and related uncertainties, and the criteria we used to classify runaway stars. In Sect. 5 we present the results we obtained for the field and runaway stars. In Sect. 6 we discuss the velocity dispersion of the field stars, the spatial distribution of the runaways, compare our results with previous works, discuss the percentage of runaways as a function of the spectral type, and briefly comment on the X-ray and gamma-ray binaries we found. We conclude by summarizing our main findings in Sect. 7.

2. Catalogs

2.1. Gaia

Gaia is an astrometric mission of the European Space Agency (ESA) that was launched in 2013. The main objective of this spacecraft is to measure the 3D spatial and velocity distribution of stars in the Galaxy with unprecedented precision. It also provides astrophysical properties of stars, including surface gravity and effective temperature. Thus, the mission allows mapping and revealing the formation, structure, and evolution of the Milky Way in the magnitude range G = 3–21 (Gaia Collaboration 2016).

Over time, Gaia data are made public through data releases that are opened to the scientific community. The latest release, Gaia DR3 (Gaia Collaboration 2023b), was made public on 13 June 2022 and is used in this work. This release, which contains more than 1.8 × 10⁹ objects, constitutes a breakthrough in terms of quality, quantity and variety of astrophysical data. In particular, it represents the largest data set of all sky spectrophotometry and radial velocities. However, although radial velocities for 33.8 × 10⁶ stars are included, these velocities are scarce and sometimes inaccurate for the massive stars of our interest because of the priors that were used (see, e.g., Drew et al. 2022). The astrometric data of DR3, namely positions (α, δ), proper motions (μ_α^*=μ_α cos δ,μ_δ), and parallaxes (ϖ), are the same as those already included in Gaia EDR3 (Gaia Collaboration 2021b), with reference epoch 2016.0.

We point out, however, that the EDR3 parallax uncertainties were underestimated, as explained in Fabricius et al. (2021) and also noted by other authors (e.g., Maíz Apellániz et al. 2021). In addition, Maíz Apellániz (2022) provides an alternative method for deriving parallax zeropoints to the methods proposed by Lindegren et al. (2021), which is better suited for bright stars. Therefore, to obtain accurate and precise parallaxes we followed the procedure presented in Maíz Apellániz et al. (2021) and Maíz Apellániz (2022). We note that this procedure requires G > 6 and five- or six-parameter solutions in Gaia. However, we also used these restrictions when we applied our own quality cuts (see Sect. 3). To apply this procedure, first we computed corrected parallaxes using the equation ϖ_c = ϖ − Z_EDR3, where ϖ are the original Gaia parallaxes and Z_EDR3 is the parallax zeropoint, which depends on magnitude, color, and ecliptic latitude, using Eq. (4) of Maíz Apellániz (2022). Second, we computed the so-called external parallax uncertainties using their formula , where k is a magnitude-dependent constant, σ_int are the original Gaia parallax uncertainties, and σ_s is the systematic uncertainty for individual parallaxes.

2.2. Galactic O-Star Catalog

The Galactic O-Star Catalog (GOSC) is an ongoing project that collects the most accurate information of Galactic O-type stars (Maíz Apellániz et al. 2013). In this catalog, special attention is paid to the correct spectral classification of the stars. To achieve this goal, this team started the Galactic O-Star Spectroscopic Survey (GOSSS; Maíz Apellániz et al. 2011). GOSSS contains blue-violet spectra of Galactic O stars with intermediate spectral resolution (R ∼ 2500) and high signal-to-noise ratio (≥300). Accordingly, from version three onwards, GOSC contains GOSSS spectral types. GOSC provides J2000 coordinates, spectral types, luminosity classes, and related clusters and associations, among many other items. The astrometric precision is 0.001 s in right ascension and 0.01″ in declination. The V magnitude ranges from ∼2 to ∼14.

In this work, we use the main catalog and supplement 2 of GOSC, whose current versions (v4.2) contain 611 O stars and 32 BA stars, respectively, with a total of 643 stars. This catalog contains precise spectral classifications for all stars, allowing us to select only OB-type stars. Since for this catalog we are only interested in early-type stars up to spectral type B1, we removed two A0 stars from the catalog. We also eliminated 37 sources that were identified with the multiple star system flag because their astrometric data might not be reliable enough (not only in GOSC, but presumably also in Gaia). Consequently, the catalog contains 604 stars.

2.3. Be Star Spectra

The Be Star Spectra (BeSS) Database (Neiner et al. 2011) provides a catalog of Be stars that is as complete as possible. It assembles spectra obtained by professional and amateur astronomers. The current version contains 2330 Be stars, Herbig Ae/Be stars, and B[e] supergiants of the Galaxy and the Large and Small Magellanic Clouds (LMC and SMC). The catalog provides J2000 coordinates of each star and the apparent magnitude V and spectral type for most of them. The precision is 0.01 s in right ascension and 0.01″ in declination. For this catalog, the visual magnitude V ranges from ∼2 to ∼20. The information is collected from the Simbad database and the literature, when this is available. Therefore, although this catalog may contain observational biases, it includes current information of a large number of Be stars, which makes it useful for our purposes. We eliminated all stars with spectral type A or that were flagged as Herbig stars from BeSS, which represented 88 and 55 stars, respectively. We also removed 126 and 215 stars with coordinates closer than 10° to the centers of the LMC and SMC, respectively (van der Marel 2001), because we are only interested in Galactic runaways. As a result, the catalog contained 1881 stars (some stars fulfilled more than one requirement for their exclusion). We note that a few Oe stars are present in both the GOSC and BeSS catalogs, but we did not remove them from any of them despite the overlap because we wished to analyze both catalogs independently.

3. Cross match of the catalogs

3.1. Cross match of GOSC and BeSS with Gaia DR3

We first cross-matched the GOSC catalog with Gaia DR3 using a radius of 1″ within the Gaia cross match tool of the ESA Gaia archive¹. We computed the differences between the GOSC and Gaia DR3 coordinates and found standard deviations of ∼007 in both coordinates. We restricted the cross-match radius to 05 to avoid false positives, which represents about seven standard deviations. Of the 604 GOSC stars, we obtained 600 cross matches.

In the case of BeSS, we proceeded in a similar way and found standard deviations of ∼017. However, the distributions looked Gaussian in the center, but showed long tails. We therefore decided to use a cross-match radius of 15 for BeSS, which represents about nine standard deviations and includes correct matches of stars in the long tails. We obtained 1863 cross matches for the 1881 BeSS stars. Because the BeSS coordinates might not be very accurate, we compared the BeSS V and Gaia DR3 G magnitudes to search for possible misidentifications. We inspected all cases where |V − G|> 1 and verified that the Gaia source identifier or source_id provided by the Simbad database was the same as the one provided by the Gaia cross-match tool. This was always the case, except for a single star, Cl* NGC 6871 BP 2, which had |V − G|> 4, for which we manually reassigned the correct Gaia source_id.

The cross matches yielded stars with more than one Gaia DR3 source_id. For the GOSC catalog, we obtained two stars with two different Gaia DR3 source_id, while for the BeSS catalog, we obtained 55 stars. Of these, we only retained the source_ids whose coordinates were closer to those of the respective catalogs, which in all cases had similar magnitudes. After the cross match, we had 598 GOSC stars (6 stars were lost) and 1808 BeSS stars (73 stars were lost).

We then applied different quality cuts to the resulting catalogs to ensure a good quality of the astrometric data. These quality cuts are summarized below and are presented in Table A.1, together with the number of stars affected by each of them. The first cut corresponds to the opposite case of the previous paragraph: the same Gaia DR3 source_id attributed to two different stars. We found 6 stars in pairs with the same source_id for the GOSC catalog and 2 stars for the BeSS catalog. We removed all of them because the accuracy of the astrometric solution would be compromised. We then applied the quality cuts recommended by Lindegren et al. (2021) for Gaia data, as well as some additional cuts used by other authors (e.g., Bailer-Jones 2015). In this way, we rejected stars without a five- and six-parameter solution, a G magnitude lower than 6, fewer than ten visibility periods, an any apparent fractional parallax error (ratio of the corrected external parallax uncertainty to the corrected parallax, σ_ext/ϖ_c) higher than 0.2, a negative parallax, and a renormalized unit weight error (RUWE) parameter > 1.35, which indicates a poor astrometric solution. We used a slightly lower value of RUWE than 1.4 that is mentioned in the Gaia data release documentation¹ to ensure good-quality data. As we explain in Sect. 4, we also excluded stars with galactocentric distances smaller than 5 kpc, which led to an additional loss of 6 GOSC stars and 1 BeSS star. After applying all these cuts, we had what we call the GOSC-Gaia DR3 catalog, which contains 417 stars, and the BeSS-Gaia DR3 catalog, which contains 1335 stars (we note that 12 of these are Oe stars that are also present in the GOSC-Gaia DR3 catalog).

3.2. Properties of GOSC-Gaia DR3 and BeSS-Gaia DR3

Figure 1 shows the distribution of G magnitudes for the stars of the GOSC-Gaia DR3 catalog in blue and the BeSS-Gaia DR3 catalog in red. The distributions are represented using a bin size of 0.5 magnitudes. The lower limit of the distribution corresponds to the G > 6 quality cut, and the upper limit reaches ∼15 and ∼17 for the O- and Be-type stars, respectively. The magnitudes of the GOSC-Gaia DR3 stars have a maximum around 8–9, and the BeSS-Gaia DR3 stars have maxima of ∼9 and ∼14. The upper end of the first maximum is due to the observational limit of the instruments used for the majority of the Be star studies. The second maximum corresponds to a detection of Be stars with larger telescopes or very low-resolution surveys (Neiner et al. 2011).

Fig. 1. Histogram of the G magnitudes for GOSC-Gaia DR3 stars in blue and BeSS-Gaia DR3 stars in red, with a bin size of 0.5 magnitudes.

To compute distances d to the stars (from the Sun), we inverted the corrected Gaia parallaxes ϖ_c discussed in Sect. 2. This procedure introduces significant biases and provides asymmetric distance uncertainties, particularly for fractional parallax uncertainties above our quality-cut limit of 0.2 (Bailer-Jones 2015). In addition, we are aware that the quality cut in fractional parallax uncertainty biases the sample toward nearby and bright objects (Luri et al. 2018). However, our original sample is already formed by bright and relatively nearby objects, and we only lost an additional ∼4% of the stars by using this quality cut in addition to the others discussed above. In contrast, this allowed us to obtain a cleaner sample. We show the distribution of the distances from the Sun for both catalogs in Fig. 2. The stars are mainly located up to 3 kpc from the Sun, but the distributions extend up to 9 kpc. The nearest GOSC-Gaia DR3 stars are at ∼0.7 kpc because of the G > 6 magnitude limit, which prevents the presence of nearby bright stars. The BeSS-Gaia DR3 stars are typically closer. Their distribution starts at ∼0.1 kpc because they are fainter and, considering the whole sample, less affected by the magnitude limit.

Fig. 2. Histogram of the distances to the Sun of GOSC-Gaia DR3 stars in blue and of BeSS-Gaia DR3 stars in red, with a bin size of 0.25 kpc.

To display the positions of the stars in our catalogs in galactocentric Cartesian coordinates, we used the convention in which the X-axis starts at the Galactic center and points toward the Sun, the Y-axis is perpendicular to it and defined positive toward the direction of the Galactic rotation for the Sun, and the Z-axis is perpendicular to the X- and Y-axes. To compute these coordinates, we first computed the (l, b) coordinates using the two equatorial coordinates of the North Galactic Pole (NGP), (α_NGP,δ_NGP), and the position angle of the North Celestial Pole with respect to the great semicircle passing through the zero Galactic longitude and the NGP, θ₀, all taken from Reid & Brunthaler (2004). Second, we used the obtained (l, b) angular coordinates and considered a galactocentric solar radius of R_⊙ = 8.15 kpc (Reid et al. 2019) to compute the galactocentric Cartesian coordinates XYZ and the galactocentric distance R_*. We show in Fig. 3 the XY galactocentric coordinates for the stars in our catalogs, where the Galactic center is located at (0,0) kpc. The absence of nearby GOSC-Gaia DR3 stars due to the cut in magnitude is clearly visible compared to the distribution of BeSS-Gaia DR3 stars. O-type stars, which are younger than Be stars, are expected to follow the spiral arms of the Milky Way better. Comparing Fig. 3 with Fig. 2 of Xu et al. (2021), we can see that our O-type stars approximately trace the Perseus, the Local, and the Sagittarius-Carina arms. In addition, a comparison of Fig. 3 with Fig. 5 of Pantaleoni González et al. (2021) reveals that O-type stars also lies within the two interarm structures they discovered: the Vela OB1 association and the Cepheus spur. The overdensity of BeSS-Gaia DR3 stars in the direction from the Sun toward XY = (13, 5) kpc is due to low-resolution surveys of Be stars in that particular direction, and it corresponds to the second magnitude peak in Fig. 1. Although not shown here, most of the stars in our catalogs are located close to the Galactic plane, with mean and standard deviations of Z = −0.00 ± 0.11 kpc for the GOSC-Gaia DR3 stars and Z = −0.01 ± 0.19 kpc for the BeSS-Gaia DR3 stars.

Fig. 3. Galactocentric XY coordinates. The coordinates of the GOSC-Gaia DR3 stars are shown in blue, and those of the BeSS-Gaia DR3 stars are plotted in red. The position of the Sun is marked with a yellow circle at (X,Y) = (8.15,0) kpc. The Galactic center is located at (0,0) kpc.

4. Search for runaways: Method

Our goal is to search for runaway stars among the stars in our catalogs. Runaway stars move with a significant peculiar velocity with respect to the mean Galactic rotation, which is the one followed by field stars. To reach this goal, we have to compute and analyze the velocities of the stars, handle the propagation of uncertainties, and define a criterion to classify the stars as runaways.

4.1. Velocities with respect to the regional standard of rest

The Galactic velocity components for any given star are defined as follows: U is positive toward the Galactic center, V is positive in the direction of the Galactic rotation, and W is positive toward the NGP. To compute the velocities of the stars, we used the different reference systems presented in Fig. 4. First, we computed the Galactic velocity components with respect to the local standard of rest (LSR), (U_LSR,V_LSR, and W_LSR). To do this, we applied the transformations of Johnson & Soderblom (1987) and added the solar motion with respect to the LSR (U_⊙,V_⊙,W_⊙). In these equations, we took the equatorial coordinates (α,δ) and the proper motions (μ_α^*=μ_α cos δ,μ_δ) from Gaia DR3, our corrected Gaia DR3 parallaxes ϖ_c (see Sect. 3), and (α_NGP,δ_NGP) and θ₀ from Reid & Brunthaler (2004). We used the solar motion values provided by the A5 fit of Reid et al. (2019): U_⊙ = 10.8 ± 1.8, V_⊙ = 13.6 ± 6.8, and W_⊙ = 7.6 ± 1.0 km s⁻¹.

Fig. 4. Reference systems used in this work, depicted in the XY plane. The Sun and the Galactic center are indicated with a yellow and black circle, respectively. The white star represents a given star in our catalogs. l is the Galactic longitude, and θ and β are defined in the text. Θ⊙ is the circular rotation velocity at the position of the Sun, and Θ*, RSR is the rotational velocity at the position of the star. LSR and RSR velocities are shown in light blue and green, respectively. The tangential velocity VTAN and the line-of-sight velocity VLOS are shown in red. The third velocity components of the LSR and the RSR systems, WLSR and WRSR, are not shown in the figure because they are perpendicular to this plane.

The transformations of Johnson & Soderblom (1987) require the radial velocity ρ, which is not available in Gaia DR3 data for most of our stars (only ∼10% of the GOSC-Gaia DR3 stars, and ∼3% of the BeSS-Gaia DR3 stars have Gaia radial velocities). To circumvent this problem, we assigned the theoretical radial velocity to each start that it should have if it were moving in its regional standard of rest (RSR) according to a Galactic rotation curve with a velocity Θ_*, RSR at galactocentric distance R_*. We used the Galactic rotation curve provided by the A5 fit of Reid et al. (2019): circular rotation velocity at the position of the Sun Θ_⊙ = 236 ± 7 km s⁻¹, and a distance from the Sun to the Galactic center of R_⊙ = 8.15 ± 0.15 kpc. Because these authors did not provide the rotation curve slope, we applied a linear fit to their data listed in Table 4, Col. 3, which are shown in the bottom panel of their Fig. 11. However, the Galactic rotation curve deviates from the linear behavior for galactocentric distances smaller than 5 kpc, therefore we discarded a few stars with R_* < 5 kpc, as already mentioned in Sect. 3. Our linear fit for galactocentric distances greater than 5 kpc, fixing Θ = Θ_⊙ for R = R_⊙, provides a slope of dΘ/dR = −0.4 ± 0.4 km s⁻¹ kpc⁻¹. With R_⊙, Θ_⊙, dΘ/dR, and R_* we computed Θ_*, RSR for each star and derived the theoretical radial velocity considering the Galactic coordinates (l, b), the angle θ of the star shown in Fig. 4, Θ_⊙ and U_⊙, V_⊙, and W_⊙. This allowed us to determine U_LSR, V_LSR, and W_LSR.

The next step was to determine the velocities of the star with respect to the RSR. The velocities were computed as follows:

(1)

The distributions of the GOSC-Gaia DR3 and BeSS-Gaia DR3 velocities (U_RSR,V_RSR,W_RSR) are presented and discussed in Appendix B.

At this point, we defined a new reference system by applying a rotation to the Galactic plane velocities U_RSR and V_RSR (see Fig. 4) in such a way that the unknown radial velocity mainly affects the line-of-sight component, V_LOS, and does not affect the other component, V_TAN. This allowed us to minimize the influence of the theoretical radial velocity. The components of the new system were computed as follows:

(2)

where V_TAN is the tangential velocity contained in the plane of the sky and parallel to the Galactic plane, V_LOS is the line-of-sight velocity, and W_RSR did not change. In these equations, β = l + θ − 90. From then on, we left V_LOS and worked with the 2D system composed of V_TAN and W_RSR.

The histograms of the GOSC-Gaia DR3 velocities V_TAN and W_RSR with a bin size of 2 km s⁻¹ are shown in Fig. 5 (the V_TAN axis is inverted according to the geometry of our reference system). The figures are limited to the ±100 km s⁻¹ abscissa range for better visualization. Two and three stars are beyond this range for the V_TAN and the W_RSR components, respectively, and their highest absolute velocities are |V_TAN| = 179.7 and |W_RSR| = 129.1 km s⁻¹. The histograms show velocities clustered around zero that approximately follow Gaussian functions (superimposed in the figures), and some stars clearly deviate from this behavior. This is expected for field stars and for runaway stars, respectively. We note that stars with V_TAN = 0 follow the Galactic rotation curve, while stars with W_RSR = 0 only move in the Galactic plane. The Gaussian fits shown in the figure were computed considering all the stars and are thus affected by the runaway stars, which are identified in the following sections. The histograms of the BeSS-Gaia DR3 velocities V_TAN and W_RSR with a bin size of 2 km s⁻¹ are shown in Fig. 6. The figures are limited to the ±100 km s⁻¹ abscissa range for better visualization. Only one star lies outside this range in the V_TAN component, whose absolute velocity is |V_TAN| = 131.1 km s⁻¹.

Fig. 5. Distributions of the VTAN and WRSR velocities for the GOSC-Gaia DR3 stars. The histograms have a bin size of 2 km s−1. The orange lines represent Gaussian functions (GF) fitted to the data, whose means and standard deviations in km s−1 are quoted above the panels. The abscissa ranges have been limited to ±100 km s−1. Top: histogram of the VTAN velocities. Bottom: histogram of the WRSR velocities.

Fig. 6. Same as Fig. 5 but for the BeSS-Gaia DR3 stars and with the Gaussian functions fitted to the data shown in blue.

4.2. Propagation of uncertainties

In all our computations, the sources of uncertainty are the five astrometric parameters (α,δ,μ_α^*,μ_δ, and ϖ_c), the solar velocities (U_⊙,V_⊙, and W_⊙), the distance from the Sun to the Galactic center R_⊙, the circular rotation velocity at the position of the Sun Θ_⊙, and the slope of the Galactic rotation curve dΘ/dR. We call X any of the parameters of the sources of uncertainty (α,δ,μ_α^*,μ_δ,ϖ_c,U_⊙,V_⊙,W_⊙,R_⊙,Θ_⊙, and dΘ/dR), and Y the functions in which the uncertainties are propagated: positions, angles, and velocities. The error on a function Y was computed as , where J_X is the Jacobian matrix of all the parameters X, and C is the covariance matrix, which only contains diagonal elements in the case of independent variables². These diagonal elements are the uncertainties of each of our parameters σ_X. This method assumes symmetric uncertainties, as is the case here.

The uncertainties in V_TAN for the GOSC-Gaia DR3 catalog range between 1.5 and 35.8 km s⁻¹, while 95% of the stars have uncertainties smaller than 7.1 km s⁻¹. For W_RSR, these numbers are 0.7, 24.9, and 4.7 km s⁻¹. The uncertainties in V_TAN for the BeSS-Gaia DR3 catalog range between 1.3 and 17.1 km s⁻¹, while 95% of the stars have uncertainties below 6.1 km s⁻¹. For W_RSR, these numbers are 0.7, 14.0, and 2.0 km s⁻¹. We show in Fig. 7 the 95th percentile of the V_TAN and W_RSR velocity uncertainties for different source of uncertainty and considering all the uncertainties together (which correspond to the numbers listed above). The different sources of uncertainty considered are the proper motions (μ_α^* and μ_δ), the corrected parallax (ϖ_c), the solar motion (U_⊙,V_⊙,W_⊙), and the uncertainties associated with the Galactic rotation curve (R_⊙,Θ_⊙,dΘ/dR). Uncertainties in position are not shown because they were negligible and had no influence on the velocity uncertainties. Figure 7 shows that the uncertainties in proper motion produce velocity uncertainties around 0.4 km s⁻¹ in all cases. The velocity uncertainties for V_TAN are dominated by the uncertainties in the solar motion (particularly due to V_⊙ = 13.6 ± 6.8 km s⁻¹) in the case of BeSS-Gaia DR3, while an additional contribution from the uncertainty in corrected parallax takes place for GOSC-Gaia DR3. The contribution from the uncertainties related to the Galactic rotation curve are only about 2–3 km s⁻¹. The velocity uncertainties for W_RSR are dominated by the uncertainties in the corrected parallax owing to the small uncertainty of the solar motion perpendicular to the Galactic plane (W_⊙ = 7.6 ± 1.0 km s⁻¹) and to the lack of influence of uncertainties of the Galactic rotation curve in the velocity perpendicular to the Galactic plane. Considering all the uncertainties together, the V_TAN uncertainties are larger than the corresponding W_RSR uncertainties mainly because of the uncertainties related to the corrected parallax, the solar motion, and the Galactic rotation curve. We also note that the GOSC-Gaia DR3 uncertainties are generally larger than those of the BeSS-Gaia DR3. The reason is that O stars lie at larger distances on average than Be stars and that the applied correction to their parallaxes and corresponding uncertainties are also larger (because of their magnitudes and colors according to Maíz Apellániz 2022).

Fig. 7. Ninety-fifth percentile of the VTAN and WRSR velocity uncertainties for different sources of uncertainty, and considering all the uncertainties together, for GOSC-Gaia DR3 stars in blue and BeSS-Gaia DR3 stars in red.

4.3. Field stars versus runaway stars

Runaway stars have historically been classified as stars that move with a significant peculiar velocity, either in two or three dimensions, that is, also above a certain threshold. This threshold was initially placed at 40 km s⁻¹ (Blaauw 1961). However, it was reduced to 30 km s⁻¹ by Gies (1987) based on a radial velocity component that is 3 times the standard deviation of the residual peculiar radial velocity for the bulk of O stars in the sample. This threshold has been further reduced to ∼25–28 km s⁻¹ by more recent works (Tetzlaff et al. 2011; Kobulnicky & Chick 2022). This reduction is based on the results of Tetzlaff et al. (2011), who found that the peculiar 2D tangential velocity had two components described by two Maxwellian distributions: a low- and high-velocity component, which intersect at 20 km s⁻¹. It is therefore sufficient to place the threshold at around ∼25 km s⁻¹ to separate the stars of the high-velocity group. In addition, Boubert & Evans (2018) found a median runaway velocity for their Be stars of 19.6 km s⁻¹, which is even lower than the ∼25 km s⁻¹ threshold. This would imply a significant loss of their runaways. Stone (1991) already noted that it would be better to search for runaways with a method independent of a runaway velocity threshold, and they suggested to do this by fitting a bimodal distribution function to the velocity histograms. Our samples are not large enough to properly trace the distribution for runaway stars. However, we used Gaussian functions to characterize the velocity distributions of field stars, which allowed us to classify the stars with significantly higher velocities as runaways without the need of a velocity threshold.

In this work, we define a runaway star as an star with a high peculiar 2D velocity (V_TAN,W_RSR) with respect to the mean Galactic rotation at a 3-sigma confidence level (but see below). This was achieved through the following parameter:

(3)

where and are the means of the velocity distributions obtained from the Gaussian fits, and and are velocity uncertainties computed using the following equations: and . These last equations take the standard deviations of the velocity distributions obtained from the Gaussian fits and into account, as well as the individual uncertainties of the velocities of each star and . We note that in all this procedure, we evaluated how significant the V_TAN and W_RSR velocities of the stars are with respect the mean velocities defined by the stars themselves (not strictly with respect to the Galactic rotation curve). The stars with E > 1 were classified as runaways. However, the means and standard deviations obtained from the Gaussian fits, particularly the latter, are affected by the presence of runaway stars. Therefore, we needed to proceed in an iterative way using a 3-sigma clipping algorithm as follows:

1. Using Eq. (3), we identified all the stars with E > 1.

2. We produced histograms without these stars and performed new Gaussian fits in both velocities.

3. Using Eq. (3) with updated values of , , , and , we identified stars with E > 1.

4. When we found more stars that fulfilled E > 1, we returned to point 2, otherwise, we stopped.

After this clipping process, the means and standard deviations of the Gaussian fits represent the distribution of the field stars better. At this point, stars with E > 1 were classified as runaways, and stars with E < 1 were classified as field stars.

5. Search for runaways: Results

In this section, we present our results for the velocity and spatial distributions of the field and the runaway stars. We also present the results of an additional search for runaways using our method for specific spectral type bins.

5.1. Velocity and spatial distributions

We show in Table 1 the means and standard deviations for the field stars after clipping of the Gaussian fits applied to the V_TAN and W_RSR velocity distributions (using a bin size of 2 km s⁻¹), and of the Z and b distributions for the complete and nearby (d < 2 kpc) GOSC-Gaia DR3 and BeSS-Gaia DR3 catalogs. We note that is larger than for both catalogs, with a higher value of for the BeSS- stars than for those from GOSC-Gaia DR3. In all cases, the mean values are closer to zero than the bin size. The mean Z and b values are close to zero with some dispersion.

After the clipping process discussed in Sect. 4.3, we obtained 106 runaway stars for the GOSC-Gaia DR3 catalog (42 newly published as such here, to our knowledge), which is equivalent to 25.4% of the catalog. For the BeSS-Gaia DR3 catalog, we obtained 69 runaway stars (47 newly published as such here), representing 5.2% of the catalog. We note that the already known runaway Oe star HD 57682 is also classified as runaway in the GOSC-Gaia DR3 catalog. We show in Table 2 the means and standard deviations as in Table 1, but for the runaway stars and without using Gaussian fits for the velocity distributions. As expected, the velocity dispersions are larger for the runaways than for the field stars, with values in the range 20–40 km s⁻¹. This agrees with the dispersions of 25–30 km s⁻¹ in the literature (Stone 1991; Tetzlaff et al. 2011). Finally, the dispersions in Z and b are larger than for the field stars (see below).

The 2D (V_TAN,W_RSR) velocity distribution of the GOSC-Gaia DR3 stars is displayed in Fig. 8. Runaway stars, which are indicated in blue, clearly have high velocities, whereas field stars, which are represented in black, have velocities about (0, 0) and are approximately contained in a 2D ellipse (if all velocity uncertainties of the individual stars were the same, it would be a real ellipse). We note that many of the black circles for the field stars overlap because the number of stars with low velocities is high, as shown in Fig. 5. The 2D velocity distribution of the BeSS-Gaia DR3 stars is displayed in Fig. 9. Runaway stars are indicated in red and field stars in black (the scale is different from Fig. 8). The observed behavior is similar as for the GOSC-Gaia DR3 stars, although in this case, the runaway stars have slightly lower velocities. We note that the standard deviations of the velocity distributions of the field stars, and , presented in Table 1 are fundamental for determining which stars are runaways. Higher values of the standard deviations imply a lower number of detected runaways. For both catalogs, is smaller than . In addition, as presented in Sect. 4.2, the V_TAN uncertainties are larger than those of W_RSR. This is due to the contribution of the Galactic rotation curve uncertainty, which affects the V_TAN component but not W_RSR. Both effects translate into more runaways that are found in the W_RSR component, as shown in Figs. 8 and 9, where the smaller axis of the field star ellipse is obtained in W_RSR.

Fig. 8. WRSR as a function of VTAN for the GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue.

Fig. 9. WRSR as a function of VTAN for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.

The galactocentric XY coordinates for the GOSC-Gaia DR3 and the BeSS-Gaia DR3 stars are shown in Figs. 10 and 11, respectively. As noted when presenting Fig. 3, there is an absence of GOSC-Gaia DR3 stars close to the Sun due to the cut in G magnitude. This cut also affects BeSS-Gaia DR3 stars, but the effect is not visible in the figures because of the many fainter objects. Therefore, there is a significant number of Be runaway stars close to the Sun. The figures do not show any privileged direction from the Sun with a significantly larger or smaller number of runaway stars. There is also no significant difference as a function of galactocentric distance.

Fig. 10. Galactocentric XY coordinates for the GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue. The position of the Sun is marked with a yellow circle at (X,Y) = (8.15,0) kpc.

Fig. 11. Same as Fig. 10 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.

Figures of the galactocentric XZ and YZ coordinates for both catalogs are presented in Appendix C. The mean and standard deviation in the Z coordinate for the field and the runaway stars are presented in Tables 1 and 2, respectively. For the GOSC-Gaia DR3 catalog, the dispersion in Z for the runaways is about three times larger than for the field stars. For the BeSS-Gaia DR3 catalog, this dispersion is about six times larger for the runaway stars. When limiting the catalog to d < 2 kpc, we typically obtained lower values for both field and runaway stars.

The projections in the (l,b) Galactic coordinates for the GOSC- and the BeSS-Gaia DR3 stars are displayed in Figs. 12 and 13, respectively. In both cases, the stars are clustered around Galactic latitude b = 0° with some dispersion. The mean and standard deviations in b for the field and the runaway stars are presented in the last column of Tables 1 and 2, respectively. For the GOSC-Gaia DR3 catalog, the dispersion in b for the runaways is about twice larger than for the field stars. For the BeSS-Gaia DR3 catalog, this dispersion is about three times larger for the runaway stars. The standard deviation for the BeSS-Gaia DR3 stars is larger than for the GOSC-Gaia DR3, indicating that the last stars are more clustered and closer to the Galactic plane. When limiting the catalog to d < 2 kpc, we obtain similar values in all cases.

Fig. 12. Galactic (l,b) coordinates for GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue.

Fig. 13. Galactic (l,b) coordinates for BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.

In Tables 3 and 4 we show ten rows of data corresponding to the runaways with the higher E values, in decreasing order, of the GOSC and BeSS-Gaia DR3 catalogs, respectively. They include their names and identifications, coordinates, distances, magnitudes, spectral types, velocities, E values, a column indicating the reference in the cases that they had already been reported as runaways, and a last column with additional information. In these tables, we include the 2D peculiar velocity, which is computed as . For the GOSC-Gaia DR3 runaway stars, it ranges from 16 to 210 km s⁻¹ and has a median (mean) of ∼40 (∼46) km s⁻¹. For the BeSS-Gaia DR3 runaway stars, ranges from 16 to 132 km s⁻¹ and has a median (mean) of ∼30 (37) km s⁻¹. The complete versions of these tables, including all classified runaways with additional parameters and precision, are available at the CDS.

Finally, we note that the bin size of the V_TAN and W_RSR histograms is relevant when we apply Eq. (3) to classify the runaway stars. In this work, we used a bin size of 2 km s⁻¹. However, we tried with a bin size of 3 km s⁻¹ and the results did not change in the classification of the GOSC-Gaia DR3 runaways. For the BeSS-Gaia DR3 catalog, we only lost two runaways. The mean and standard deviations obtained from the Gaussian fits did not change significantly, except for of the BeSS-Gaia DR3 catalog, which changed from 9.3 to 9.5 km s⁻¹.

5.2. Runaways as a function of spectral type

We performed an additional search for runaways using our method for specific spectral type bins. We first removed stars from the GOSC-Gaia DR3 and BeSS-Gaia DR3 catalogs with an uncertain or variable spectral classification and binary systems with two spectral classifications. Next, we divided our samples into different spectral type bins. We tried to choose the bins as homogeneously as possible, considering that the stars are not homogeneously distributed in spectral type and that the GOSC and BeSS samples are different in size. We decided to use only two bins for each catalog to obtain significantly large samples (plus an additional bin in the intersection). To obtain the number of runaway stars, we applied Eq. (3) for each bin independently. The results in number and percentage are shown in Table 5. These results are similar to those obtained with the full samples although there is a trend of a decreasing percentage of runaway stars at later spectral types, as shown in Fig. 14. Because the percentage drops very significantly from the O to the Be stars, we conducted the search for an intermediate bin of spectral type O8–B1e composed of about 200 GOSC and 200 BeSS stars. The obtained results are shown in the last row of Table 5. They reveal a percentage between those found previously.

Fig. 14. Percentage of runaway stars as a function of spectral type. Blue and red correspond to the GOSC-Gaia DR3 and BeSS-Gaia DR3 stars, respectively.

6. Discussion

This section is organized into five subsections. First, we discuss the velocity dispersion of the field stars. Second, we examine the spatial distribution of the runaway stars. Third, we compare our results with previous works and discuss the implications of our findings. Fourth, we investigate the percentage of runaways as a function of spectral type and its relation to the origin of the runaway stars. Finally, we place our findings in the context of high-mass X-ray binaries and gamma-ray binaries.

6.1. Velocity dispersion of field stars

The kinematics of OB field stars is related to that of the gas of the Milky Way disk. Because these stars are young, they should approximately follow the motion of the gas from which they recently formed (Gontcharov 2012; Gaia Collaboration 2023a). In this context, we compared our velocity dispersions with those of two works. On the one hand, Gaia Collaboration (2023a) presented an average azimuthal and vertical dispersions of 9.4 and 6.2 km s⁻¹, respectively, for a sample of OB stars. Our dispersions of ∼7–9 km s⁻¹ for V_TAN and ∼5 km s⁻¹ for W_RSR (see Table 1) are slightly smaller. This may be because here we present the dispersions of a clean sample of stars without runaway stars, which would widen the distributions. We note that although the azimuthal velocity is not the same as our tangential velocity, their dispersions are expected to be comparable (see Eq. (2)). On the other hand, Marasco et al. (2017) obtained a velocity dispersion of the gas of about 4–9 km s⁻¹, which agrees very well with both Gaia Collaboration (2023a) and our results.

As shown in Table 1, the dispersion of W_RSR is around 5 km s⁻¹ for both catalogs, while the dispersion in V_TAN increases to 6.6 km s⁻¹ for GOSC-Gaia DR3 stars and up to 9.3 km s⁻¹ for BeSS-Gaia DR3 stars. Therefore, there is a velocity dispersion increase in the Galactic plane. This can be due to a combined effect of parallax (distance) uncertainties and the application of the Galactic rotation curve, which induces errors for distant stars in the tangential component of the velocity, but not in the vertical component. To study this possibility, we performed the same study, but restricted the distance to 2 kpc. The results, also shown in Table 1, reveal that the dispersion of V_TAN for GOSC-Gaia DR3 stars decreases to 5.0 km s⁻¹ and is thus very similar to the different dispersions obtained in W_RSR. It is therefore clear that the use of the Galactic rotation curve for stars with larger distances (and with larger uncertainties) artificially enhances this dispersion. However, in the case of BeSS-Gaia DR3 stars, the dispersion of V_TAN does not decrease and seems to be intrinsic. This can be understood because Be stars are older than O stars and are therefore more affected by Galactic velocity diffusion in the disk, producing a larger velocity dispersion (see, e.g., Fig. 10 of Gaia Collaboration (2021a) for trends in the same sense). We conclude that the velocity dispersion at birth for young stars is ∼5 km s⁻¹, while it increases to higher values with age due to Galactic velocity diffusion. Earlier studies claiming velocities of ∼10 km s⁻¹ for young stars were probably affected by not so accurate measurements and their influence in the computation of velocities (e.g., Mihalas & Binney 1981 or Tetzlaff et al. 2011).

6.2. Spatial distribution of the runaway stars

The galactocentric XY coordinates of the runaway stars of the GOSC- and BeSS-Gaia DR3 catalogs were shown in Figs. 10 and 11, respectively. As noted in Sect. 5.1, there is no preferred spatial distribution for the runaways as there is for field stars, but there are no O-type runaway stars close to the Sun, while Be runaways are clustered around the solar neighborhood. To explain this behavior, two factors are named. On the one hand, because the brightest stars with G < 6 were filtered, nearly all O-type stars were removed around the solar neighborhood, while this was not the case for the Be stars. On the other hand, the closer stars have smaller uncertainties and are therefore easier to be classified as runaways if they have a significant velocity. These are the reasons for the many runaway Be stars close to the Sun, while this is not the case for O stars.

The distribution of runaway stars might not trace the spiral arms because massive runaway stars are moving with high velocities and have been expelled from their birthplace (Blaauw 1993; Maíz Apellániz et al. 2018). Therefore, these stars can be found outside the spiral arms and are in particular more spread out along the Z-axis or the Galactic latitude b than the field stars. This is shown for the Z-axis in Figs. C.1 and C.2. Field stars are mostly located close to the Galactic plane, while runaway stars are more distributed along the Z-axis. A similar effect is shown for the Galactic latitude b in Figs. 12 and 13. The dispersions we obtained for Z and b for runaways (see Table 2) are significantly larger than for field stars (see Table 1). For the GOSC-Gaia DR3 catalog, the dispersion in Z for the runaways is about three times larger than for the field stars. This factor is about six for the BeSS-Gaia DR3 catalog. For b, these numbers are about two and about three, respectively. The dispersion differences in Z and b for the O and Be runaways are larger for Be runaways. This is explained by the fact that the Be stars are older, which allowed them to wander farther from the Galactic plane.

We note that the mean and standard deviation in Z and b were also obtained for stars up to 2 kpc only. For the field stars (Table 1), there are no significant differences with the complete catalog results. This indicates that our results are stable and that for the BeSS-Gaia DR3 catalog, they are not biased by the overdensity toward XY = (13, 5) kpc that extends from Z ≃ 0 kpc close to the Sun to Z ≃ 0.4 kpc away from it. In the case of the runaway stars (see Table 2), the dispersions in b do not change significantly, while there is a small decrease in the dispersions in Z. This last issue arises because runaway stars with higher values of |Z| are rare and also located at larger XY distances from the Sun, which effectively removed them with our d < 2 kpc filter.

6.3. Comparison with previous works

For each of our runaways, we searched in the literature whether they had already been identified as runaways by other authors. As presented in Sect. 5.1, we found 42 and 47 runaway stars with no previous identifications for the GOSC- and the BeSS-Gaia DR3 catalogs, respectively. We recovered 64 and 22 known runaway stars for the two catalogs, respectively. The corresponding publication references are listed in the column “Run Ref.” of Tables 3 and 4. Here we mention the most common publications for our runaways. Considering both catalogs, we recovered 15, 41, 23, and 12 runaway stars that were previously identified by Tetzlaff et al. (2011), Maíz Apellániz et al. (2018), Kobulnicky & Chick (2022), and Boubert & Evans (2018), respectively. We note that Kobulnicky & Chick (2022) only published the names of 25 out of the 102 GOSC runaways they identified. Therefore, we cannot compare our results with their complete list. We show a summary of the samples, methods, thresholds and percentage of runaways published in these works in Table 6, together with the results obtained in this work.

Tetzlaff et al. (2011) searched for young (≤50 Myr) runaways around 3 kpc from the Sun in three dimensions using HIPPARCOS data and radial velocities from the literature. They found a runaway percentage of ∼27%. The limit in age used by these authors includes basically our first three bins of Table 5, for which we would obtain ∼13.5% of runaways, but using a 2D method. Although the samples are not the same and might be biased, we seem to recover only roughly 50% of the runaway stars because we did not use radial velocities. Therefore, the percentages shown in Table 5 should probably be multiplied by a factor ∼2 to obtain the real percentages of runaways, which would imply percentages above 50% for O-type stars. Because this percentage is quite high, we performed simulations of runaway stars using a 3D velocity method that reproduced our 2D results. We obtained a 3D percentage of ∼30% for O-type stars compared to the 2D percentage of 25.4% (see Appendix D).

Our GOSC-Gaia DR3 runaway percentage of 25.4% (see Sect. 5.1) is significantly higher than that of Maíz Apellániz et al. (2018), who obtained 5.7% for O stars. Differences lie slightly in the samples, but mainly in the Gaia data, DR1 by them versus DR3 here, and especially in the method used, 2D proper motions by them versus 2D space velocities here. This prevented them from identifying stars with high peculiar velocity at large distances.

Kobulnicky & Chick (2022) searched for runaway stars within GOSC using Gaia EDR3 data and a 2D method, following the procedure of computing the velocities from Randall et al. (2015). There are some notable differences between their approach and ours. Kobulnicky & Chick (2022) did not use the zeropoint correction to the parallax, and most notably, the larger external parallax uncertainties proposed by Maíz Apellániz (2022). Therefore, their results could be too optimistic in this sense. We also used a different rotation curve and solar motion values to compute the velocities. Although they used a Monte Carlo analysis to estimate the final uncertainties, they considered as input uncertainties those in parallax, proper motions, and solar motion, while we also considered the uncertainty in the galactocentric distance of the Sun and in the Galactic rotation curve (plus the small coordinate uncertainties). Finally, they classified as runaway stars those with a peculiar 2D velocity above a threshold 25 km s⁻¹. Of the 25 GOSC runaways with published names, we recover only 23, mainly because of our larger but more realistic external parallax uncertainties. We also note that for these stars, we obtain slightly lower velocities than theirs on average because of the parallax correction we applied (the use of different Galactic rotation curves and solar motion velocities had practically no influence). The total number and percentage of runaways obtained by these authors are very similar to ours: 102 by them versus 106 here, that is, ∼22% versus ∼25%. However, the final samples are different: They considered as runaway stars those with higher velocities because of their 25 km s⁻¹ threshold, although their uncertainties are underestimated, while we used more realistic uncertainties, but had a threshold of 16 km s⁻¹ that was derived from the dispersion of the field stars and the individual uncertainties. This means that our sample contains fewer false positives related to uncertainties, but includes stars that were classified as walkaways by other authors.

The above works are mainly concerned with runaway O stars. Although studies of B runaway stars can be found in the literature (Hoogerwerf et al. 2001; Tetzlaff et al. 2011), those dealing specifically with Be stars are less common. Boubert & Evans (2018) exclusively searched for Be runaway stars within the BeSS and LAMOST (Hou et al. 2016) catalogs and the stars of Berger & Gies (2001) that fulfill several quality cuts, using Gaia DR1 data and radial velocities from the literature (3D method). To classify the stars as runaways, they used a Bayesian method in which they compared the hypothesis that every star in their catalog is a runaway with the null hypothesis that it is a member of the Milky Way disk using priors. Within their sample of 632 Be stars, they found 13.1% runaway stars, representing 83 objects. We recovered 12 of them in our BeSS-Gaia DR3 runaways. Their percentage of runaways of 13.1% is considerably higher than the one obtained here of 5.2% and those found in the literature for B stars (see Sect. 6.4). This discrepancy of a factor ∼2.5 could be understood considering that they used a 3D Bayesian method and Gaia DR1 astrometric data. Again, we performed simulations of runaway stars using a 3D velocity method that reproduced our 2D results and obtained a 3D percentage of ∼6.7% for Be stars compared to the 2D percentage of 5.2% (see Appendix D). Berger & Gies (2001) also analyzed Be runaway stars with a 3D method and found ∼6.7% runaway stars.

6.4. Runaways as a function of the spectral type

The percentage of runaways as a function of the spectral type that was presented in Sect. 5.2 decreases as follows: 25.1%, 23.7%, 6.2%, and 4.8% for O2−O7, O8−O9, B0e−B3e, and B4e−B9e, respectively. This trend was already noted by Blaauw (1961), who showed from 3D velocities that the runaway frequency as a function of the spectral type drops sharply from roughly 20% among O-type stars to 2.5% among B0–B0.5 stars, and it is even lower, 1.5%, among B1–B5 stars. Gies & Bolton (1986) obtained lower limits and estimated similar percentages of OB runaway stars, while Stone (1979) increased the percentage of runaway O stars to ∼49%. More recently, Tetzlaff et al. (2011) found that ∼27% of the young stars are runaways, and Boubert & Evans (2018) provided 13.1% runaways in a sample of Be stars. The scatter among the results is therefore considerable. It mainly falls on the different input data, the method (2D or 3D), and the existence or absence of a velocity threshold to classify the stars as runaways. However, there is agreement in the sense that the percentage of runaway O stars is significantly higher than for B or Be stars.

In this context, two possible scenarios might explain the origin of runaway stars: the dynamical ejection scenario (DES; Poveda et al. 1967), and the binary supernova scenario (BSS; Blaauw 1961). Dorigo Jones et al. (2020) compared these two scenarios for OB runaway stars in the SMC. They reported that the DES scenario often results in faster, more massive runaways, while the BSS scenario results in lower-velocity ejections of lower-mass stars. Their results for the SMC show that the dynamical mechanism, which favors massive runaways, dominates by a factor of 2–3. The GOSC-Gaia DR3 stars have higher velocities in general than those in BeSS-Gaia DR3 (see Figs. 8 and 9), which agrees with the scenarios discussed above. Moreover, we obtain a factor of ∼5 between the percentage of runaway stars among O-type stars versus Be-type stars in the Galaxy. This reinforces the dominance of the DES scenario versus the BSS one discussed by Dorigo Jones et al. (2020).

6.5. High-mass X-ray and gamma-ray binaries

We cross-checked our catalogs of runaway stars with the catalog of High Mass X-ray Binaries (HMXBs) by Fortin et al. (2023), which also contains gamma-ray binaries. For the GOSC-Gaia DR3 runaways, we found the following five objects in common with our names (their names) ordered by decreasing values of the E parameter: V479 Sct (LS 5039), GP Vel (Vela X-1), HD 153919 (4U 1700−377), LM Vel (IGR J08408−4503), and Cygnus X-1. For the BeSS-Gaia DR3 runaways, we found the following two objects in common: GSC 03588-00834 (SAX J2103.5+4545) and BQ Cam (V 0332+53). LM Vel and GSC 03588-00834 are identified here as runaway HMXBs for the first time.

Of all these HMXBs, we would like to comment on the particular case of Cygnus X-1. We obtained a value of E = 1.10, which is close to our threshold criterion to classify stars as runaways. The study by Mirabel & Rodrigues (2003) revealed that this object is moving with a space velocity of only 9 ± 2 km s⁻¹ with respect to Cyg OB3. They used this information to justify that Cygnus X-1 is a member of Cyg OB3 association and has a low space velocity with respect to its environment. However, the stars of Cyg OB3 are already moving with a space velocity of ∼22 km s⁻¹ with respect to the mean Galactic rotation curve and basically toward the Galactic center (Rao et al. 2020), which in this case nearly coincides with the V_TAN direction. Therefore, the slight additional velocity of Cygnus X-1 is enough for us to classify it as a runaway object, although it would not be classified as such considering the space velocity of the association. This is a limitation of our method, which could affect some stars with values of E ≳ 1.

On the other hand, we are also interested in searching for new gamma-ray binaries. Only ten such systems are currently known according to Bordas (2023). Four of the eight systems with good spectral type classification for the massive stars are O stars and four are Be stars. Only the O star LS 5039 is listed in the current version of the GOSC catalog, and we classify it as a runaway, as was already known (Ribó et al. 2002; Moldón et al. 2012). Three Be gamma-ray binaries are listed in the current version of the BeSS catalog: PSR B1259−63, LS I +61 303, and HESS J0632+057. Of these, only PSR B1259−63 is a runaway, but mainly in the radial direction, in which we are blind. Therefore, none of them was classified as a runaway by our method, as expected. The detection of LS 5039 reveals that a careful inspection of our runaway stars, particularly among O-type stars, using multiwavelength catalogs could reveal new gamma-ray binaries.

7. Summary and conclusions

Using Gaia DR3 data, we have searched for Galactic massive runaway stars of O and Be spectral type. We cross-matched the GOSC and the BeSS catalog with Gaia DR3 and obtained 417 and 1335 stars, respectively. We presented a 2D method to classify the stars as runaways if they had a significant velocity at a 3-σ confidence level with respect to the field stars, which follow the Galactic rotation curve. The main conclusions of this work are summarized below.

• For O-type field stars, we obtained a velocity dispersion perpendicular to the Galactic plane of ∼5 km s⁻¹ and a parallel dispersion of ∼7 km s⁻¹. This increase is due to parallax uncertainties and the use of the Galactic rotation curve for distant stars.

• For Be-type field stars, we obtained a velocity dispersion perpendicular to the Galactic plane of ∼5 km s⁻¹ and a parallel dispersion of ∼9 km s⁻¹. This increase is dominated by the Galactic velocity diffusion allowed by their older ages.

• We found 106 O runaway stars, representing 25.4% of our O-type star catalog. Forty-two of them were not previously identified as runaways.

• We found 69 Be runaway stars, representing 5.2% of our Be-type star catalog. Forty-seven of them were not previously identified as runaways.

• Considering O- and Be-type stars, the dispersion of runaways is about three to six and about two to three times higher in Z and b than that of field stars, respectively. This is explained by the ejections they underwent when they became runaways.

• The percentage of runaways decreases drastically from O- to Be-type stars, but it decreases gradually within each type: ∼25%, ∼24%, ∼6%, and ∼5% for O2−O7, O8−O9, B0e−B3e, and B4e−B9e, respectively.

• Comparison with previous works using 3D methods suggests that these percentages could be ∼2 and ∼2.5 times higher for O- and Be-type stars, respectively. However, 3D simulations for our catalogs reveal that our 2D percentages could increase only up to ∼30% and ∼6.7%, respectively, implying factors of ∼1.2 and ∼1.3, respectively.

• We find higher velocities for O-type runaways than for Be-type runaways. This agrees with previous works. We find a factor of ∼5 between the percentage of runaway stars among O-type stars versus Be-type stars. Both facts underline that the DES scenario is more likely than the BSS scenario.

• Seven of our runaway stars are HMXBs, two of which are identified as new runaways in this work, and one is a gamma-ray binary. This opens the door to identifying new high-energy systems among our runaways by conducting detailed multiwavelength searches.

The upcoming Gaia data releases, containing significantly improved parallaxes, together with the use of radial velocities, could lead us to discover more runaway stars among our input samples.

Acknowledgments

We thank the anonymous referee for useful suggestions and comments that helped to improve the content of the manuscript. We would like to thank the members of the Gaia team at Universitat de Barcelona for many useful discussions. We thank J. Maíz Apellániz for valuable discussions. We acknowledge financial support from the State Agency for Research of the Spanish Ministry of Science and Innovation under grants PID2019-105510GB-C31/AEI/10.13039/501100011033, PID2019-104114RB-C33/AEI/10.13039/501100011033, PID2022-136828NB-C41/AEI/10.13039/501100011033/ERDF/EU, and PID2022-138172NB-C43/AEI/10.13039/501100011033/ERDF/EU, and through the Unit of Excellence María de Maeztu 2020-2023 award to the Institute of Cosmos Sciences (CEX2019-000918-M) 2021SGR00679. M.C.-C. acknowledges the grant PRE2020-094140 funded by MCIN/AEI/10.13039/501100011033 and FSE/ESF funds. 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 work has made use of the BeSS database, operated at LESIA, Observatoire de Meudon, France: http://basebe.obspm.fr. This research has made use of NASA’s Astrophysics Data System. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

References

Appendix A: Quality cuts applied to the cross matches

Table A.1 shows the quality cuts we applied to the 598 and 1808 stars obtained after the cross matches of GOSC and BeSS with Gaia DR3, respectively. Column 1 lists the quality cuts we used. Columns 2–3 show the number of GOSC stars (#) affected by each quality cut and the corresponding percentage (%) with respect to the total of 598 stars. Columns 4–5 show the same for BeSS stars, where percentages are given with respect to the total of 1808 stars. The last row of the table includes the resulting numbers after all the cuts were applied. This last row does not show the sum of the star numbers and percentages because some stars are affected by more than one quality cut. There are 423 GOSC stars and 1336 BeSS stars that survive after these quality cuts.

Appendix B: U_RSR, V_RSR, and W_RSR velocity distributions

The computed values of U_RSR, V_RSR, and W_RSR are clustered around 0 km s⁻¹. We produced histograms with a bin size of 2 km s⁻¹ to better display the velocity distributions (as a compromise between resolution and significance). The histograms of the RSR velocities of the GOSC-Gaia DR3 stars are shown in Fig. B.1. The top, middle, and bottom panels correspond to the U_RSR, V_RSR, and W_RSR velocities, respectively. The figures are limited to the ±100 km s⁻¹ abscissa range for better visualization. For the U_RSR and W_RSR velocities, two and three stars lie outside that range, respectively. The maximum absolute velocities of these three stars are |U_RSR| = 177.1 km s⁻¹ and |W_RSR| = 129.1 km s⁻¹. No stars lie outside the V_RSR velocity range. The histograms show velocities clustered around zero, with some dispersion that approximately follows a Gaussian function. Some stars display much higher velocities. Gaussian fits are also shown in Fig. B.1. The same histograms are shown for the BeSS-Gaia DR3 stars in Fig. B.2. The figures are also limited to the ±100 km s⁻¹ abscissa range for better visualization. Only one star lies outside that range in U_RSR, with |U_RSR| = 129.6 km s⁻¹, while no stars lie outside this range for the V_RSR and W_RSR velocities.

Fig. B.1. Distributions of the URSR, VRSR, and WRSR velocities for the GOSC-Gaia DR3 stars. The histograms have a bin size of 2 km s−1. The orange lines represent Gaussian functions fitted to the data, whose means and standard deviations in km s−1 are quoted above the panels. The abscissa ranges have been limited to ±100 km s−1. Top: Histogram of the URSR velocities. Middle: Histogram of the VRSR velocities. Bottom: Histogram of the WRSR velocities.

Fig. B.2. Same as Fig. B.1, but for the BeSS-Gaia DR3 stars and with the Gaussian functions fitted to the data shown in blue.

For the GOSC- and BeSS-Gaia DR3 catalogs, the means and standard deviations of the Gaussian fits applied to the RSR velocity distributions using a bin size of 2 km s⁻¹ are shown in Table B.1. We note that the mean values are always closer to 0 than the bin size of 2 km s⁻¹. The dispersions obtained in U_RSR, V_RSR, and W_RSR are about 5 km s⁻¹, except for V_RSR for the GOSC-Gaia DR3 stars, which is significantly smaller. We emphasize that all velocities, but particularly U_RSR and V_RSR, are affected by the lack of observational radial velocities. Our theoretical radial velocities decrease the dispersions of the distributions. For comparison, Tetzlaff et al. (2011) used proper motions from HIPPARCOS and radial velocities from the literature for 7663 young stars to obtain dispersions in U_RSR, V_RSR, and W_RSR of 10.7, 10.7, and 5.3 km s⁻¹, respectively. The dispersion in W_RSR we derived is similar to theirs. Our dispersions in U_RSR and V_RSR are significantly smaller than theirs, as expected due to the use of theoretical radial velocities.

Appendix C: Galactocentric XZ and YZ planes

Figure C.1 shows the galactocentric XZ (top) and YZ (bottom) coordinates for the stars of the GOSC-Gaia DR3 catalog. Figure C.2 shows the galactocentric XZ (top) and YZ (bottom) coordinates for the stars of the BeSS-Gaia DR3 catalog. In these figures, the Z coordinate is limited to between ±1 kpc for a better visualization. Five stars in the BeSS-Gaia DR3 catalog lie outside that range, and their highest |Z| value is 3.7 kpc. The GOSC-Gaia DR3 stars are clustered around the Galactic plane with some dispersion (see main text). The BeSS-Gaia DR3 stars show a similar behavior, although the overdensity discussed in Sect. 3 toward (X, Y) = (13, 5) is located above the Galactic plane.

Fig. C.1. Galactocentric XZ and YZ coordinates for the stars of the GOSC-Gaia DR3 catalog. Field stars are depicted in black, and runaway stars are shown in blue. The position of the Sun is marked with a yellow circle. The ordinate axis is limited to ±1 kpc and the axes have different scales. Top: Galactocentric XZ coordinates. Bottom: Galactocentric YZ coordinates, but with a different scale in the abscissa axis.

Fig. C.2. Same as Fig. C.1 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.

Appendix D: Simulations of runaway stars using a 3D velocity method

In this work, we used a 2D (V_TAN,W_RSR) velocity method to compute the number of runaway stars in our catalogs. Here we report on the estimation of the number of runaway stars that we would detect using a 3D velocity method with simulated velocity distributions that reproduce our 2D results for runaway stars. We used the following approach.

1. We used an exponential decay function to simulate the modules of the 3D peculiar velocities for a population of N stars. Specifically, we chose the python random.exponential function of the NumPy library, which draws samples from an exponential distribution following the equation f(x, β) = (1/β)exp( − x/β) with a given β parameter. In addition, we multiplied the output values by a parameter v₀, in km s⁻¹, which scales the exponential function to our desired range of velocities. We have three free parameters: the number of simulated stars N, the β parameter of the exponential function, and the velocity scale factor v₀.

2. We assumed a 3D velocity distribution with spherical symmetry. To do this, we used a uniformly distributed random angle ϕ ∈ [0,2π) and a uniformly distributed random function cos θ ∈ [−1,1], where θ ∈ [0,π].

3. We computed the projections of the 3D simulated velocities using the velocity modules and the ϕ and θ angles.

4. For the velocity uncertainties of the stars, we assumed those corresponding to the 68% percentile of the 2D uncertainty distributions of (V_TAN,W_RSR). For the component, we considered the same uncertainties as for V_TAN because both are in the Galactic plane and are affected by similar uncertainties (corrected parallax, solar motion, and Galactic rotation curve). The velocity uncertainties we assumed are presented in the last three columns of Table D.1.

5. We considered an underlying population of field stars. To do this, we took the standard deviations of (V_TAN,W_RSR) obtained from the Gaussian fits, but centered at 0 km s⁻¹ for simplicity (the differences are negligible). Again, for , we considered the same standard deviation as for V_TAN. The means and standard deviations are presented in Table D.1.

6. We classified 2D runaway stars as those with E > 1 in Eq. (3) using instead of (V_TAN,W_RSR) and the input parameters presented in Table D.1. Similarly, we also classified 3D runaway stars adding the corresponding term, , and input parameters to Eq. (3).

7. We conducted 10000 simulations to compute the average distributions of the simulated velocities , the mean numbers of 2D and 3D runaway stars, and the percentages with respect to the number of stars in our catalogs.

8. We performed the following comparisons between the simulated population of 2D runaway stars found using and the population of 2D runaway stars found using (V_TAN,W_RSR):

   • The number of runaway stars in 2D.

   • The standard deviations of the 2D velocities: , (to be compared with σ_VTAN and σ_WRSR of Table 2).

   • The histograms of , , and with respect to those of V_TAN, W_RSR, and .

   • The 2D velocity plot vs. with respect to W_RSR vs. V_TAN.

9. After these comparisons, we modified the three free parameters of the simulations in step 1 to reproduce the 2D runaway results obtained in this work as closely as possible, and we repeated the whole process until convergence was achieved. For the GOSC-Gaia DR3 catalog, convergence was achieved for the following parameters: N = 247, β = 0.14, v₀ = 230 km s⁻¹. In the case of BeSS-Gaia DR3 catalog, the values were N = 328, β = 0.11, and v₀ = 130 km s⁻¹.

The results of our simulations for O-type stars are shown in Table D.2. In the 2D case, we accurately reproduce the number and percentage of runaway stars and the standard deviations of their velocities (see Table 2). We show in the top panel of Fig. D.1 the 2D velocity distribution of the simulated runaway stars in blue for one of the 10000 simulations. The field stars of the GOSC-Gaia DR3 catalog are plotted in black for reference. This figure resembles Fig. 8. In the 3D case, we obtained 125 runaway stars, representing 30% of the 417 stars in the GOSC-Gaia DR3 catalog. The standard deviations of the runaway stars preserve the spherical symmetry, and are smaller than those in 2D because we can now detect stars as runaways with lower and velocities. This is clearly visible in the bottom panel of Fig. D.1, which includes the simulated 3D runaway stars. In this figure, runaway stars lie in the region in which field stars were located in 2D because of the contribution of the (for both positive and negative values of ).

Fig. D.1. as a function of for the simulated runaway stars represented in blue for one of the 10000 simulations that match the GOSC-Gaia DR3 results. The (VTAN,WRSR) velocities of the field stars of the GOSC-Gaia DR3 catalog are depicted in black. Top: Simulated 2D runaway stars. Bottom: Simulated 3D runaway stars with either positive or negative values of .

For the Be-type stars, we applied the same process as described in steps 1 to 9 with the input parameters presented in Table D.1, but with a slight difference. In step 3, we multiplied the and distributions by a factor 1.5 to reproduce the BeSS-Gaia DR3 velocity distributions because the standard deviation of V_TAN obtained from the Gaussian fit is clearly higher than that of W_RSR (see Table 1). We note that this can be understood as a result of Galactic velocity diffusion in the disk (see Sect. 6.1). The results of our simulations for Be-type stars are shown in Table D.2. In the 2D case, the number, percentage, and the velocity standard deviations of the runaway stars were also reproduced (see Table 2). We show in the top panel of Fig. D.2 the 2D velocity distribution of the simulated runaway stars in red for one of the 10000 simulations. The field stars of the BeSS-Gaia DR3 catalog are shown in black for reference. This figure is similar to Fig. 9. In the 3D case, we obtained 90 runaway stars, representing 6.7% of the 1335 stars in the BeSS-Gaia DR3 catalog. As expected due to the 1.5 factors used, the standard deviations of the runaway stars are equal for the and components, and smaller for the component. They are also smaller than those in 2D. The 2D velocity distribution including the simulated runaway stars in 3D is shown in the bottom panel of Fig. D.2, where we also have new runaways in the area of the 2D field stars because of the contribution.

Fig. D.2. Same as Fig. D.1 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.

When the 2D and 3D results are compared, the number of runaways for O-type stars increased by 19 ± 4 from 2D to 3D (the uncertainty reflects the standard deviation of the increase in the 10000 simulations). This represents an increase from ∼25% in 2D to ∼30% in 3D, which corresponds to a factor of ∼1.2. For the Be stars, the increase in the number of runaways was 21 ± 4 from 2D to 3D. This represents an increase from ∼5.2% in 2D to ∼6.7% in 3D, corresponding to a factor of ∼1.3.

Finally, we note that we also used a power-law function to simulate the modules of the 3D peculiar velocities in step 1 of the simulations presented above and obtained similar results, although the exponential decay function reproduced the observational results better.

All Tables

All Figures

	Fig. 1. Histogram of the G magnitudes for GOSC-Gaia DR3 stars in blue and BeSS-Gaia DR3 stars in red, with a bin size of 0.5 magnitudes.
In the text

	Fig. 2. Histogram of the distances to the Sun of GOSC-Gaia DR3 stars in blue and of BeSS-Gaia DR3 stars in red, with a bin size of 0.25 kpc.
In the text

	Fig. 3. Galactocentric XY coordinates. The coordinates of the GOSC-Gaia DR3 stars are shown in blue, and those of the BeSS-Gaia DR3 stars are plotted in red. The position of the Sun is marked with a yellow circle at (X,Y) = (8.15,0) kpc. The Galactic center is located at (0,0) kpc.
In the text

Fig. 4. Reference systems used in this work, depicted in the XY plane. The Sun and the Galactic center are indicated with a yellow and black circle, respectively. The white star represents a given star in our catalogs. l is the Galactic longitude, and θ and β are defined in the text. Θ⊙ is the circular rotation velocity at the position of the Sun, and Θ*, RSR is the rotational velocity at the position of the star. LSR and RSR velocities are shown in light blue and green, respectively. The tangential velocity VTAN and the line-of-sight velocity VLOS are shown in red. The third velocity components of the LSR and the RSR systems, WLSR and WRSR, are not shown in the figure because they are perpendicular to this plane.

In the text

Fig. 5. Distributions of the VTAN and WRSR velocities for the GOSC-Gaia DR3 stars. The histograms have a bin size of 2 km s−1. The orange lines represent Gaussian functions (GF) fitted to the data, whose means and standard deviations in km s−1 are quoted above the panels. The abscissa ranges have been limited to ±100 km s−1. Top: histogram of the VTAN velocities. Bottom: histogram of the WRSR velocities.

In the text

	Fig. 6. Same as Fig. 5 but for the BeSS-Gaia DR3 stars and with the Gaussian functions fitted to the data shown in blue.
In the text

	Fig. 7. Ninety-fifth percentile of the VTAN and WRSR velocity uncertainties for different sources of uncertainty, and considering all the uncertainties together, for GOSC-Gaia DR3 stars in blue and BeSS-Gaia DR3 stars in red.
In the text

	Fig. 8. WRSR as a function of VTAN for the GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue.
In the text

	Fig. 9. WRSR as a function of VTAN for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.
In the text

	Fig. 10. Galactocentric XY coordinates for the GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue. The position of the Sun is marked with a yellow circle at (X,Y) = (8.15,0) kpc.
In the text

	Fig. 11. Same as Fig. 10 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.
In the text

	Fig. 12. Galactic (l,b) coordinates for GOSC-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in blue.
In the text

	Fig. 13. Galactic (l,b) coordinates for BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.
In the text

	Fig. 14. Percentage of runaway stars as a function of spectral type. Blue and red correspond to the GOSC-Gaia DR3 and BeSS-Gaia DR3 stars, respectively.
In the text

Fig. B.1. Distributions of the URSR, VRSR, and WRSR velocities for the GOSC-Gaia DR3 stars. The histograms have a bin size of 2 km s−1. The orange lines represent Gaussian functions fitted to the data, whose means and standard deviations in km s−1 are quoted above the panels. The abscissa ranges have been limited to ±100 km s−1. Top: Histogram of the URSR velocities. Middle: Histogram of the VRSR velocities. Bottom: Histogram of the WRSR velocities.

In the text

	Fig. B.2. Same as Fig. B.1, but for the BeSS-Gaia DR3 stars and with the Gaussian functions fitted to the data shown in blue.
In the text

Fig. C.1. Galactocentric XZ and YZ coordinates for the stars of the GOSC-Gaia DR3 catalog. Field stars are depicted in black, and runaway stars are shown in blue. The position of the Sun is marked with a yellow circle. The ordinate axis is limited to ±1 kpc and the axes have different scales. Top: Galactocentric XZ coordinates. Bottom: Galactocentric YZ coordinates, but with a different scale in the abscissa axis.

In the text

	Fig. C.2. Same as Fig. C.1 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.
In the text

	Fig. D.1. as a function of for the simulated runaway stars represented in blue for one of the 10000 simulations that match the GOSC-Gaia DR3 results. The (VTAN,WRSR) velocities of the field stars of the GOSC-Gaia DR3 catalog are depicted in black. Top: Simulated 2D runaway stars. Bottom: Simulated 3D runaway stars with either positive or negative values of .
In the text

	Fig. D.2. Same as Fig. D.1 for the BeSS-Gaia DR3 stars. Field stars are depicted in black, and runaway stars are shown in red.
In the text