Galactic runaway O and Be stars found using Gaia DR3

A relevant fraction of massive stars are runaway stars. These stars move with a significant peculiar velocity with respect to their environment. We aim to discover and characterize the population of massive and early-type runaway stars in the GOSC and BeSS catalogs using Gaia DR3 astrometric data. We present a 2-dimensional method in the velocity space to discover runaway stars as those that deviate significantly from the velocity distribution of field stars, which are considered to follow the Galactic rotation curve. We found 106 O runaway stars, 42 of which were not previously identified as runaways. We found 69 Be runaway stars, 47 of which were not previously identified as runaways. The dispersion of runaway stars is a few times higher in Z and b than that of field stars. This is explained by the ejections they underwent when they became runaways. The percentage of runaways is 25.4% for O-type stars, and it is 5.2% for Be-type stars. In addition, we conducted simulations in 3 dimensions for our catalogs. They revealed that these percentages could increase to ~30% and ~6.7%, respectively. Our runaway stars include seven X-ray binaries and one gamma-ray binary. Moreover, we obtain velocity dispersions of ~5 km/s perpendicular to the Galactic plane for O- and Be-type field stars. These values increase in the Galactic plane to ~7 km/s for O-type stars due to uncertainties and to ~9 km/s for Be-type stars due to Galactic velocity diffusion. The excellent Gaia DR3 astrometric data have allowed us to identify a significant number of O-type and Be-type runaways in the GOSC and BeSS catalogs. The higher percentages and higher velocities found for O-type compared to Be-type runaways underline that the dynamical ejection scenario is more likely than the binary supernova scenario. Our results open the door to identifying new high-energy systems among our runaways by conducting detailed studies.


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 7 yr), while B stars have longer life-spans (up to ∼10 8 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), gammaray binaries (Dubus 2013), millisecond pulsar systems, and dou-⋆ The complete versions of Tables 3 and 4 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr(130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/vol/page ⋆⋆ Serra Húnter Fellow ble 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 Way in the magnitude range G = 3-21 (Gaia Collaboration et al. 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 et al. 2023b), was made public on 13 June 2022 and is used in this work.This release, which contains more than 1.8 × 10 9 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 6 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 et al. 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) andMaí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.

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.

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.

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 1 .We computed the differences between the GOSC and Gaia DR3 coordinates and found standard deviations of ∼0 ′′ .07 in both coordinates.We restricted the cross-match radius to 0 ′′ .5 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 ∼0 ′′ .17.However, the distributions looked Gaussian in the center, but showed long tails.We therefore decided to use a cross-match radius of 1 ′′ .5 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 re-1 gea.esac.esa.int/archivespective 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 1 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).

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).
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.
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, θ 0 , 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 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.

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 θ 0 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 −1 .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 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 V TAN and the line-of-sight velocity V LOS are shown in red.The third velocity components of the LSR and the RSR systems, W LSR and W RSR , are not shown in the figure because they are perpendicular to this plane.

∼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 −1 , 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, Column 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 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: 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: 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-ofsight 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 −1 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 −1 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 −1 .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 −1 are shown in Fig. 6.The figures are limited to the ±100 km s −1 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 −1 .

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 J X CJ T X , 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 variables2 .These diagonal elements are the uncertainties of each of our pa- rameters σ 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 −1 , while 95% of the stars have uncertainties smaller than 7.1 km s −1 .For W RSR , these numbers are 0.7, 24.9, and 4.7 km s −1 .The uncertainties in V TAN for the BeSS-Gaia DR3 catalog range between 1.3 and 17.1 km s −1 , while 95% of the stars have uncertainties below 6.1 km s −1 .For W RSR , these numbers are 0.7, 14.0, and 2.0 km s −1 .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 −1 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 −1 ) 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 −1 .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 −1 ) 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).

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 −1 (Blaauw 1961).However, it was reduced to 30 km s −1 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 −1 by more recent works (Tetzlaff et al. 2011 ) nicky & 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 −1 .It is therefore sufficient to place the threshold at around ∼25 km s −1 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 −1 , which is even lower than the ∼25 km s −1 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: where µ GF V TAN and µ GF W RSR are the means of the velocity distributions obtained from the Gaussian fits, and σV TAN and σW RSR are velocity uncertainties computed using the following equations: W RSR, * .These last equations take the standard deviations of the velocity distributions obtained from the Gaussian fits σ GF V TAN and σ GF W RSR into account, as well as the individual uncertainties of the velocities of each star σ V TAN, * and σ W RSR, * .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 µ GF V TAN , µ GF W RSR , σV TAN , and σW RSR , 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.

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.

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 −1 ), 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 σ GF V TAN is larger than σ GF W RSR for both catalogs, with a higher value of σ GF V TAN 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 −1 .This agrees with the dispersions of 25-30 km s −1 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

Catalog
Stars Field Stars   slightly lower velocities.We note that the standard deviations of the velocity distributions of the field stars, σ GF V TAN and σ GF W RSR , 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, σ GF W RSR is smaller than σ GF V TAN .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 .
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.
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.
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 −1 and has a median (mean) of ∼40 (∼46) km s −1 .For the BeSS-Gaia DR3 runaway stars, V 2D PEC ranges from 16 to 132 km s −1 and has a median (mean) of ∼30 (37) km s −1 .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 −1 .However, we tried with a bin size of 3 km s −1 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 σ GF V TAN of the BeSS-Gaia DR3 catalog, which changed from 9.3 to 9.5 km s −1 .

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.

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.

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 et al. 2023a).In this context, we compared our velocity dispersions with those of two works.On the one hand, Gaia Collaboration et al. (2023a) presented an average azimuthal and vertical dispersions of 9.4 and 6.2 km s −1 , respectively, for a sample of OB stars.Our dispersions of ∼7-9 km s −1 for V TAN and ∼5 km s −1 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 −1 , which agrees very well with both Gaia Collaboration et al. (2023a) and our results.
As shown in Table 1, the dispersion of W RSR is around 5 km s −1 for both catalogs, while the dispersion in V TAN increases to 6.6 km s −1 for GOSC-Gaia DR3 stars and up to 9.3 km s −1 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 −1 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 et al. (2021a) for trends in the same sense).We conclude that the velocity dispersion at birth for young stars is ∼5 km s −1 , while it increases to higher values with age due to Galactic velocity diffusion.Earlier studies claiming velocities of ∼10 km s −1 for young stars were probably affected by not so accurate measurements and their influence in the computation of velocities (e.g., Mihalas & Binney 1981or Tetzlaff et al. 2011).

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 Sec.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  Kobulnicky & Chick (2022).Comments column contains additional information: 'f' for fast rotator stars (Britavskiy et al. 2023); 'b' for binary hints (Brandt 2021), 'h' for high-mass X-ray binaries (Fortin et al. 2023); 'g' for gamma-ray binaries (Bordas 2023).The complete version of this table, including the 106 classified runaways with additional parameters and precision, is available at the CDS .
Table 4. Data of the 10 BeSS-Gaia DR3 runaway stars with the higher values of E in decreasing order.
BeSS Name  & Evans (2018).Comments column contains additional information: 'f' for fast rotator stars (Britavskiy et al. 2023); 'b' for binary hints (Brandt 2021), 'h' for high-mass X-ray binaries (Fortin et al. 2023); 'g' for gamma-ray binaries (Bordas 2023).The complete version of this table, including the 69 classified runaways with additional parameters and precision, is available at the CDS.
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.

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   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 Apel-lá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 stars with a peculiar 2D velocity above a threshold 25 km s −1 .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 −1 threshold, although their uncertainties are underestimated, while we used more realistic uncertainties, but had a threshold of 16 km s −1 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.

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 percent-  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).

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 −1 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 −1 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.

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-sigma 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 −1 and a parallel dispersion of ∼7 km s −1 .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 −1 and a parallel dispersion of ∼9 km s −1 .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 Betype 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 Oto 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 Betype 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 gammaray binary.This opens the door to identifying new highenergy 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.
Table B.1.Means and standard deviations of the Gaussian fits applied to the RSR velocity distributions using a bin size of 2 km s −1 .

Catalog
Stars

Fig. 1 .Fig. 2 .
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.

Fig. 3 .
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.

Fig. 4 .
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 V TAN and the line-of-sight velocity V LOS are shown in red.The third velocity components of the LSR and the RSR systems, W LSR and W RSR , are not shown in the figure because they are perpendicular to this plane.

Fig. 5 .
Fig. 5. Distributions of the V TAN and W RSR 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 V TAN velocities.Bottom: histogram of the W RSR velocities.

Fig. 6 .
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.

Fig. 7 .
Fig. 7. Ninety-fifth percentile of the V TAN and W RSR 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.

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

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

Fig. 10 .Fig. 11 .
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. 14 .
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.
Fig. B.1.Distributions of the U RSR , V RSR , and W RSR 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 U RSR velocities.Middle: Histogram of the V RSR velocities.Bottom: Histogram of the W RSR velocities.

Fig. C. 1 .
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 .
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.

Table 1 .
Means and standard deviations of different distributions for the field stars after clipping.

Table 2 .
Means and standard deviations of different distributions for the runaway stars after clipping.

Table 3 .
Data of the 10 GOSC-Gaia DR3 runaway stars with the higher values of E in decreasing order.

Table 5 .
Stars and runaway stars as a function of spectral type.

Table 6 .
Blaauw 1961l.1967)18)fworks on all-sky searches of runaways among young stars discussed in Sec.6.3.More recently,Tetzlaff et al. (2011)found that ∼27% of the young stars are runaways, andBoubert & 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).DorigoJones et al. ( Table D.1.Input parameters for the simulations.
Table D.2.Output results of the simulations.