© ESO 2019

1. Introduction

Hypervelocity celestial objects have characteristic velocities of the order of 1000 km s⁻¹, meaning those travelling 100 kpc in less than 100 Myr, and they may be ubiquitous in the Universe. The existence of hypervelocity stars was first predicted by Hills (1988) and later discussed by Yu & Tremaine (2003), but the first bona fide hypervelocity star (HV 1 or SDSS J090744.99+024506.9 with a heliocentric radial velocity of 853 ± 12 km s⁻¹) was identified as such by Brown et al. (2005). The first hypervelocity star cluster, HVGC-1 or H70848 with an offset from the systemic velocity of the Virgo Cluster >2300 km s⁻¹, was found by Caldwell et al. (2014). In addition, the picture emerging from the spectroscopic analysis of large samples of star-forming galaxies in the local Universe hints at a significant number of high-velocity runaway stars and hypervelocity stars that may be following radial trajectories, away from their original hosts (for example, see Cicone et al. 2016).

It is however not easy to separate true hypervelocity stars from the high-velocity tail associated with the stellar halo, both in the Milky Way galaxy (Deason et al. 2019) and elsewhere. In the following, we will consider that hypervelocity stars are those that are not gravitationally bound to the Milky Way galaxy or to any other galaxy for that matter, though some researchers use the term to refer exclusively to those stars that were ejected after interacting with the supermassive black hole in the Galactic centre (see the review by Brown 2015). An extensive and timely updated source of data on this subject is the Open Fast Stars Catalog¹ (Guillochon et al. 2017; Boubert et al. 2018).

The known hypervelocity stars with a probable origin in the Milky Way galaxy exhibit a statistically significant anisotropic spatial spread (Brown et al. 2009) and it is not well understood why. Five recent studies (Boubert et al. 2018; Bromley et al. 2018; Hattori et al. 2018; Irrgang et al. 2018; Marchetti et al. 2019) use Gaia Data Release 2 (DR2) data (Gaia Collaboration 2016, 2018a) to independently revisit the topic of hypervelocity stars in the Milky Way. Both Boubert et al. (2018) and Irrgang et al. (2018) pay particular attention to reanalysing previously known candidates using the new data; meaning that they do not identify new ones in the Gaia catalogue.

Boubert et al. (2018) conclude that although known early-type hypervelocity stars may have been ejected from the Galactic centre, the vast majority of late-type candidates could be currently bound to the Galaxy and therefore could be part of the high-velocity tail associated with the stellar disc and halo, or perhaps could be debris from disrupted galaxies. Bromley et al. (2018) single out a number of promising new candidate hypervelocity stars at distances in the range of 10–15 kpc. Hattori et al. (2018) focus on identifying nearby, metal-poor, and relatively old extreme velocity candidates. Irrgang et al. (2018) find that among the known hypervelocity stars, the fraction with an origin in the disc rather than the Galactic centre dominates. Marchetti et al. (2019) isolate a sample that includes both candidates of Galactic and extragalactic origin, but they do not find any candidate consistent with an origin in the Galactic centre. However, they find strong evidence for a population of high-velocity stars with a probable origin outside the Milky Way galaxy. This interpretation agrees well with the ideas put forward by Cicone et al. (2016) for example. The purportedly fastest star in the Gaia catalogue, Gaia DR2 5932173855446728064 (Bromley et al. 2018; Marchetti et al. 2019), has been found to be spurious (Boubert et al. 2019).

Most new hypervelocity star candidates identified in Gaia DR2 are inside or relatively close to the nominal edge of the Milky Way disc, 15 kpc from the Galactic centre. However, stars located well beyond such a distance may have a higher probability of being hypervelocity stars and/or having an extragalactic origin. With this working hypothesis in mind, here we used Gaia DR2 data to perform a systematic exploration of the spatial distribution of distant, nominal Galactocentric distances greater than 30 kpc, hypervelocity star candidates. This investigation is aimed at understanding the origin of any putative anisotropic spatial distribution, but also at confirming or refuting the presence of a significant population of high-velocity stars of intergalactic provenance. This paper is organized as follows. Section 2 discusses data selection issues and the overall approach applied in our study. The results of the full sample are presented and discussed in Sect. 3. Section 4 focusses on the subsample with the lowest uncertainties. In Sect. 5, we evaluate the statistical significance of our findings. Results are discussed in Sect. 6 and conclusions are summarized in Sect. 7.

2. Data selection and processing

The basic premise that guides our exploration using Gaia DR2 data is that the farther the candidate hypervelocity stars are, the more likely they are to be unbound from the Milky Way galaxy; this is particularly true for candidates probably located well beyond the Galactic centre. This assumption leads us to deal with samples that have been neglected by previously published studies (Boubert et al. 2018; Bromley et al. 2018; Hattori et al. 2018; Irrgang et al. 2018; Marchetti et al. 2019) because the values of their relative parallax errors (σ_π/π, where π is the measured value of the astrometric parallax and σ_π, its uncertainty) are inherently large. It is however possible to obtain statistically significant results if the effects of the parallax errors are well characterised and understood when the sample is magnitude-limited (for example, see the discussion in Sect. 3.6 of Binney & Merrifield 1998). In the following, averages, standard deviations, medians, interquartile ranges (IQRs), and other statistical parameters have been computed in the usual way (for instance, see Press et al. 2007; Wall & Jenkins 2012).

2.1. The sample and the uncertainties

In order to interpret our results with confidence (by using colour-magnitude diagrams for example), we focussed on those sources with estimated values of the line-of-sight extinction A_G and reddening E(G_BP − G_RP); Gaia DR2 includes 87 733 672 such sources², all of them have strictly positive values of the parallax. Out of this sample, 4 831 731 sources have positions, parallax, radial velocity, and proper motions. This smaller sample is suitable for an analysis in the Galactocentric frame of reference.

We computed Galactocentric positions using the value of the distance between the Sun and the Galactic centre (Sgr A^*) given by GRAVITY Collaboration (2019), 8.18 kpc. These positions are in the Galactocentric standard of rest that is a right-handed coordinate system centred at the Galactic centre with positive axes in the directions of the Galactic centre, Galactic rotation, and the North Galactic Pole (NGP) as discussed by Johnson & Soderblom (1987) for example. Galactocentric Galactic velocity components were calculated as described by Johnson & Soderblom (1987), considering the values of the Solar motion computed by Schönrich et al. (2010) and the value of the in-plane circular motion of the local standard of rest around the Galactic centre discussed by Reid et al. (2014). These values are also used by Bromley et al. (2018). From the smaller sample, our software pipeline produced 15 681 sources with nominal Galactocentric distances >30 kpc and full kinematics with uncertainties; this is our primary full sample.

As the input data have large uncertainties, we did not use the error expressions presented by Johnson & Soderblom (1987) but a Monte Carlo (Metropolis & Ulam 1949; Press et al. 2007) sampling technique to perform error estimation. The Monte Carlo methodology used here also provides estimates for the most probable values of the various parameters. For each source in our primary sample, we generated 10⁵ realizations of its astrometric and photometric parameters to produce lists of ordered data values. Neither mean values nor standard deviations of computed parameters (with the exception of those of input Gaia DR2 data) are used in this work. As measures of central tendency and dispersion, we used the median (50th percentile) and the 16th and 84th percentiles, respectively. These are the values given in the sections, figures, and tables. The covariance matrix of parallax and proper motions was computed using the mean values, standard deviations, and correlation coefficients provided by Gaia DR2.

For example, if r is a vector made of univariate Gaussian random numbers (components r_i with i = 1, 3), the required multivariate Gaussian random samples are given by the expressions:

(1)

where π, μ_α, and μ_δ are the values of absolute stellar parallax and proper motions in right ascension and declination directions provided by Gaia DR2 and the a_ij coefficients are given by:

(2)

where σ_π, σ_μα, and σ_μδ are the standard errors in parallax and proper motions from Gaia DR2, and ρ_πμα, ρ_πμδ, and ρ_μαμδ their respective correlation coefficients, also from Gaia DR2. For each π_c, we computed the value of the distance by applying the usual relationship, d_c = 1/π_c. If the value π_c drawn from the normal distribution turned out to be negative, it was not discarded. The resulting list of 10⁵ values of d_c was ordered, and used to obtain the median and the 16th and 84th percentiles (see Fig. 1). An equivalent approach was used to estimate the relevant values of the percentiles of interest for other parameters such as the Galactocentric distance and Galactic velocity components; the values are quoted in terms of the median of the distribution, with uncertainties derived from the 16th and 84th percentiles. Given the fact that negative values of π_c were not discarded, the uncertainties were fully characterised because the information on the tail of the probability distribution of the parallax below zero was not lost in the Monte Carlo sampling.

Fig. 1. Analytical versus Monte Carlo distances for 100 representative synthetic data with σπ/π ∈ [0.2, 4.2] and π ∈ [0.0125, 0.0333] mas (uniformly distributed). These are shown with (left) and without (right) error bars, computed as described in the text. The discontinuous line shows the diagonal (d = dc) for reference.

In order to generate random numbers (r_i) with a standard Gaussian or normal distribution with mean 0 and standard deviation 1, we applied the Box–Muller method (Box & Muller 1958; Press et al. 2007). The Box–Muller method requires a uniform random variable as seed; when a computer is used to produce a uniform random variable, it will inevitably have some inaccuracies because there is a lower boundary on how close numbers can be to 0. For a 64 bits computer the smallest non-zero number is 2⁻⁶⁴, which means that the Box–Muller method will not produce random variables more than 9.42 standard deviations from the mean (see Press et al. 2007).

2.2. Parallax zero-point

Lindegren et al. (2018) use quasars to study the large-scale systematics in the astrometric quality of Gaia DR2 data, concluding that there is a significantly negative global zero-point, −0.029 mas. In other words, the reported Gaia DR2 parallaxes are too small and hence the distances are systematically overestimated if one adopts the usual relationship. In our particular case, this issue might lead to overestimated distances and therefore unrealistically high values of the Galactocentric velocity. A possible solution to this problem could be in using d_c = 1/(π_c + 0.029 mas) for the calculations described above. This is however not advisable for a number of reasons.

On the one hand, the actual value of the global zero-point is controversial. Although Graczyk et al. (2019) find a best value of the parallax zero-point of −0.031 ± 0.011 mas, which is fully consistent with the one reported by Lindegren et al. (2018), other authors find larger parallax offsets. Riess et al. (2018) give a value of −0.046 ± 0.013 mas, Leung & Bovy (2019) find a value of −0.052 ± 0.002 mas, Zinn et al. (2019) obtain −0.053 ± 0.003 mas, Schönrich et al. (2019) find a value of −0.054 ± 0.006 mas, Arenou et al. (2018) use star clusters to obtain −0.067 ± 0.012 mas (see their Fig. 16), Xu et al. (2019) compute a parallax offset of −0.075 ± 0.029 mas, and Stassun & Torres (2018) find a value of −0.082 ± 0.033 mas. The wide range in the values proposed for this offset mainly arises from the different types of objects (and thus different colour and magnitude intervals) used to compute the various estimates. Although all the proposed values of the global zero-point are below 0.100 mas in absolute terms (which is the pre-launch, expected global systematics), it must be emphasized that they are consistently negative and below the one determined by Lindegren et al. (2018), −0.029 mas. Our sample has nominal parallaxes <0.033 mas; if we perform this correction, even if we use the most optimistic value, −0.029 mas, the size of the resulting sample would become virtually zero (see Appendix A).

On the other hand, Arenou et al. (2018) show that the parallax offset is partly dependent on the scanning pattern of the Gaia mission, which makes it a function of the coordinates (see their Fig. 15). Both Lindegren et al. (2018) and Arenou et al. (2018) explicitly discourage a global zero-point correction, particularly in cases (like ours) where the sample is not well distributed over the entire sky (see below). For these reasons, we do not correct the parallaxes for the global zero-point.

2.3. Selection criteria

Our identification of distant hypervelocity star candidates was based on the values of Galactocentric distance and velocity. At Galactocentric distances >30 kpc, a Galactocentric velocity >500 km s⁻¹ is probably signalling the presence of robust hypervelocity star candidates or at least members of the high-velocity tail associated with the stellar halo – Deason et al. (2019) find a value for the escape speed in the neighbourhood of the Sun of km s⁻¹, therefore the value at distances greater than 30 kpc from the Galactic centre is expected to be lower. These criteria are more conservative than the ones used recently to discuss two new hypervelocity stars from the LAMOST spectroscopic surveys (Huang et al. 2017).

As for the origin of the candidates, we used the angle between the position and velocity vectors as a flag (this technique is discussed by Bromley et al. 2018). A hypervelocity star coming from the region of the Galactic centre must have a value of the angle close to 0° as both vectors are nearly parallel, bound stars or those moving tangentially are expected to have values of the angle in the neighbourhood of 90° as both vectors are nearly perpendicular, and those stars coming straight from intergalactic space and headed for the region of the Galactic centre should have a value of the angle close to 180° when both vectors are almost antiparallel, but in general >90°. However, the angle criterion has some weaknesses. A fast-moving star coming from outside the Milky Way galaxy may have experienced a hyperbolic encounter with the Galactic centre and the angle between the position and velocity vectors could be close to 0° when moving away from the Galactic centre and back into intergalactic space. On the other hand, any intergalactic star traversing the Galaxy without experiencing significant deflection will mimic an origin within the disc if located well outside the Galactic bulge. However, even in such cases, the estimated age of the candidate may help in discarding (or not) a possible intergalactic provenance.

The sample of high-velocity stars may be contaminated by stars with poor astrometric solutions. In order to single out those sources with reliable astrometry, we applied the criteria discussed by Lindegren et al. (2018) and Marchetti et al. (2019) (regarding the parameters astrometric_gof_al, astrometric_excess_noise_sig, mean_varpi_factor_al, visibility_periods_used, and rv_nb_transits). The data and the criteria are presented in Fig. 2. After the application of these criteria, the primary sample of 15 681 was reduced to 393 stars, 174 of them have median Galactocentric distance greater than 30 kpc and median Galactocentric Galactic velocity greater than 500 km s⁻¹ (computed using the Monte Carlo methodology described above).

Fig. 2. Quality summary of the primary sample (15 681 sources). Position and relative parallax error as a function of the parameter astrometric_excess_noise_sig or significance of excess noise (left panels). Lindegren et al. (2018) indicate that astrometrically well-behaved sources have a significance of excess noise ≤2 (region on the left of the red line). The green line in the left-hand, bottom panel signals σπ/π = 0.2. Dependency between quality flags (right panels). Lindegren et al. (2018) state that sources with astrometric_gof_al or goodness of fit <3 (below the green line of the top panel), mean_varpi_factor_al or mean parallax factor AL in the interval [−0.23, 0.32] (between the green lines of the middle panel), and visibility_periods_used or visibility periods >8 (above the green line of the bottom panel) should be used. Not shown here, but also applied, is rv_nb_transits or number of observations used to compute the value of the radial velocity >5.

2.4. Extinction and reddening

Our Gaia DR2 sources have estimated values of the line-of-sight extinction A_G and reddening E(G_BP − G_RP). These values are highly biased due to the model grids that were used to train the machine learning algorithms employed and the extinction was computed considering the usual relationship between parallax and distance, without taking into account the parallax error or parallax zero-point offset (Gaia Collaboration 2018a). Andrae et al. (2018) indicate that the typical accuracy in A_G is of the order of 0.46 mag and that of E(G_BP − G_RP) could be 0.23 mag. Classical sources to obtain estimates of Galactic dust reddening and extinction are Schlegel et al. (1998) and Schlafly & Finkbeiner (2011). Their datasets are available from the NASA/IPAC InfraRed Science Archive³ and we retrieved the relevant values to compare with those given by Gaia DR2 (although the wavelengths are different).

Figure 3 is a comparison between the extinction (left-hand side panels) and reddening (right-hand side panels) values from Gaia DR2 and those derived by Schlegel et al. (1998), top panels, and Schlafly & Finkbeiner (2011), bottom panels, for the 15 681 sources in our primary sample (in grey) and the 393 stars with the best astrometric solutions (in black). The values derived by Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) are linearly correlated, with those from Schlegel et al. (1998) being somewhat higher as seen from the reference red diagonal line, meaning E(G_BP − G_RP) = E(B − V). The error bars (only for 393 sources) of Gaia DR2 extinction and reddening show the 16th and 84th percentiles, respectively, those of Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) display the standard deviations, which are very small for reddening and zero for extinction.

Fig. 3. Extinction and reddening, a comparison. Left panels: extinction values AV from Schlegel et al. (1998), (top panels), and Schlafly & Finkbeiner (2011), (bottom panels), versus the line-of-sight extinction AG from Gaia DR2 for the 15 681 sources in our primary sample (in grey) and the 393 stars with the best astrometric solutions (in black). Right panels: colour excesses E(B − V) from Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) versus E(GBP − GRP) from Gaia DR2 for the same sources. The reference line in red is the diagonal. See the text for additional details.

Although the wavelengths are different and the uncertainties in the Gaia DR2 values are large, the agreement is fairly good for the bulk of sources with the best astrometric solutions although Gaia DR2 extinction and reddening values tend to be larger than those from Schlegel et al. (1998) and Schlafly & Finkbeiner (2011). The two-dimensional maps of Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) were obtained without considering any distance information, while the Gaia DR2 extinction values do take into account the distance derived from the usual relationship between parallax and distance, and neglecting its uncertainty. We did not observe, among the 393 stars with the best astrometric solutions, significant numbers of outliers with abnormally low or high values of extinction and reddening; therefore, it seems unlikely that our conclusions could be affected by erroneous extinction and reddening values from Gaia DR2.

2.5. Are they really that far away?

In principle, the central problem with our analysis is in the use of sources with large relative parallax errors, larger than those discussed in the five recent studies that use Gaia DR2 data (Boubert et al. 2018; Bromley et al. 2018; Hattori et al. 2018; Irrgang et al. 2018; Marchetti et al. 2019). This choice, which is the result of the basic premise pointed out above, has two main side effects: distances could be both incorrect and underestimated, and their uncertainties cannot be properly computed.

On the one hand, it is a well-known fact that when σ_π/π ≥ 0.2 it becomes very difficult to estimate the true value of the astrometric parallax (see Trumpler & Weaver 1953; Bailer-Jones 2015); in other words, the value of the distance obtained by applying the usual relationship, d = 1/π, as well as its associated error, Δd = σ_π/π², become very unreliable. Within the context of Gaia DR2 data, our self-imposed lower limit for the nominal Galactocentric distance automatically implies that our primary full sample has σ_π/π >  0.2. On the other hand, the Eddington-Trumpler-Weaver bias (Eddington 1913, 1940; Trumpler & Weaver 1953) affects parallax-limited samples such as ours. The circumstances that lead to the Eddington-Trumpler-Weaver bias have been revisited multiple times (see Lutz & Kelker 1973; Oudmaijer et al. 1998; Smith 2003; Francis 2014) and it is more commonly discussed under the term Lutz-Kelker bias; it is still unclear whether or not it can be corrected (and if it can, how, see Lutz & Kelker 1973; Koen 1992; Francis 2014). It is however widely accepted that it causes measured parallaxes to be too large, leading to underestimated distances, when the uncertainty in the measured value becomes a significant fraction of the measurement itself. In our case, this implies that the actual distance to the source is very probably larger than the one computed using the conventional expression, d = 1/π.

The analysis presented by Trumpler & Weaver (1953) suggested that any values of the distance estimated using d = 1/π for a magnitude-limited sample with σ_π/π >  0.2 could be systematically underestimated, but it did not provide a procedure to obtain the true values. At this point, one may argue that the correct approach to estimate the true distances should involve the application of Bayesian inference, adopting a prior distribution for distances in our Galaxy (see Bailer-Jones 2015; Astraatmadja & Bailer-Jones 2016; Bailer-Jones et al. 2018; Luri et al. 2018). This is particularly important for the stars in our primary full sample, since they are thought to be distant (>30 kpc), faint stars with large uncertainties on parallax. Bailer-Jones et al. (2018) compute the distances to 1.33 billion stars in Gaia DR2 using Bayesian inference and all the sources in our primary full sample have estimated distances in this catalogue. However, the fraction of distant stars in Bailer-Jones et al. (2018) is very negligible, with 31 298 sources with Heliocentric distance >15 kpc, 49 with distance >25 kpc, and 4 with distance >30 kpc; in other words, most sources with values of the parallax under 0.033 mas in Gaia DR2 have Bayesian distance estimates in Bailer-Jones et al. (2018) well below 15 kpc. The distance estimates in Bailer-Jones et al. (2018) are biased to lower values if the parallax data are poor because they adopt a small scale length for their prior. This fact interpreted within the context of the Eddington-Trumpler-Weaver bias discussed above led us to enquire how reliable were the estimates provided by Bailer-Jones et al. (2018) and (which is far more important) how strong was the effect of the Eddington-Trumpler-Weaver bias on our sample.

Eclipsing binaries are routinely used to find the distances to galaxies both in the Local Group and beyond as they provide an independent and precise alternative to astrometric parallaxes (for example, see Bell et al. 1991). In order to measure the strength of the Eddington-Trumpler-Weaver or Lutz-Kelker bias on our sample, we retrieved results from recent eclipsing binary surveys in the Magellanic Clouds. Using eclipsing binaries, Pietrzyński et al. (2013) found a value for the distance to the Large Magellanic Cloud (LMC) of 49.97 ± 0.19 (statistical) ±1.11 (systematic) kpc (or a parallax of 0.0208 ± 0.0005 mas); also observing eclipsing binaries, Hilditch et al. (2005) found a distance of 60.6 ± 1.0 (statistical) ±2.8 (systematic) kpc (or a parallax of 0.0165 ± 0.0010 mas) to the Small Magellanic Cloud (SMC). This distance interval covered reasonably well the range of interest for our primary full sample, which is 30–60 kpc for candidates with relative parallax error <1 (see Sect. 4). These values (of d and its associated π) can be used as the true ones for all practical purposes as their uncertainties are far lower than those from astrometric parallaxes. In addition, Muraveva et al. (2014) provide an extensive catalogue of 1768 eclipsing binaries in the LMC and the one in Pawlak et al. (2013) include 6138 eclipsing binaries in the SMC. The vast majority of these stars have counterparts in Gaia DR2 and, therefore, measured parallaxes. These values can be used to evaluate both the reliability of the relevant results in Bailer-Jones et al. (2018) within the context of our study and the effect of the Eddington-Trumpler-Weaver bias on our sample.

Figure 4 shows the relative differences in the values of the parallax as a function of their relative errors for sources in Muraveva et al. (2014) (top panel, LMC) and Pawlak et al. (2013) (bottom panel, SMC) with strictly positive values of the parallax in Gaia DR2. The search for matching sources was carried out using the tools provided by VizieR (Ochsenbein et al. 2000) with a radius of . The x-axis shows the value of σ_π/π from Gaia DR2, while the y-axis shows the difference between measured and true parallax divided by the uncertainty in the true parallax (meaning that for an eclipsing binary in the LMC, top panel in Fig. 4, 0.0208 mas was subtracted from the value of the Gaia DR2 parallax and divided by 0.0005 mas). Positive y-values correspond to sources with measured parallax value larger than the true one (or empirical distance lower than the real one). For sources in the LMC with strictly positive values of the parallax in Gaia DR2, the probability of measuring a value of the parallax greater than the true value is 0.83 (in other words, the probability of having a false negative if the null hypothesis is that the source is at or beyond the LMC), for the SMC we found an equivalent probability of 0.94. Within a frequentist framework and making a simple extrapolation, our interpretation of these results is that the Eddington-Trumpler-Weaver bias is present in our sample, it is very important, and over 80% of our sources may have significantly underestimated values of their distances when the usual expression, d = 1/π, is applied to convert Gaia DR2 parallaxes into distances. This bias may be even more important for sources located farther away (compare top and bottom panels in Fig. 4).

Fig. 4. Eddington-Trumpler-Weaver bias. Evaluation of this bias for data from Muraveva et al. (2014), (top panel), and Pawlak et al. (2013), (bottom panel). The y-axis shows the difference between the value of the parallax from Gaia DR2 and the one derived from the assumed distance – 0.0208 ± 0.0005 mas for the LMC (top panel, Pietrzyński et al. 2013) and 0.0165 ± 0.0010 mas for the SMC (bottom panel, Hilditch et al. 2005)— divided by the quoted uncertainty (0.0005 mas, top panel, and 0.0010 mas, bottom panel). The x-axis shows the corresponding value of the relative error in parallax from Gaia DR2. Only sources with strictly positive values of the parallax are plotted; a sample of the full catalogues is shown as the full ranges in x and y are wider. The condition σπ/π = 0.2 is plotted as a discontinuous line. Left panels: results without considering any global zero-point correction; right panels: include the correction (π + 0.029 mas) discussed in the text.

On the other hand, we repeated the analysis for several other dwarf galaxies (closer than the LMC) and obtained consistent results (virtually 100% false negatives, albeit using samples of eclipsing binaries made of just a handful of sources). If the global zero-point correction discussed above – π + 0.029 mas instead of π, meaning that some sources originally with π <  0 got positive values – is used, the effect of the Eddington-Trumpler-Weaver bias increases only slightly (see Fig. 4, right-hand side panels). The probability of having a false negative if the null hypothesis is that the source is at or beyond the LMC was 0.827 without the correction and it becomes 0.839 with the correction; for the SMC the corresponding probability goes from 0.939 to 0.943. The papers describing the samples of eclipsing binaries used in our analysis state that the fraction of foreground sources, meaning false positives, in their catalogues is negligible; therefore, the values of the probability can be considered as sufficiently reliable. In this analysis we applied a global zero-point correction to a relatively small region of the sky, although both Lindegren et al. (2018) – see their Fig. 13 – and Arenou et al. (2018) explicitly discourage this practice.

In order to evaluate the reliability of the relevant results in Bailer-Jones et al. (2018), we matched sources as pointed out above and constructed histograms with statistically meaningful bin sizes; the bin width was computed using the Freedman-Diaconis rule (Freedman & Diaconis 1981): 2 IQRn^−1/3, where n is the number of sources. Figure 5 shows the distributions of Heliocentric distances from Bailer-Jones et al. (2018) and the true distances (as vertical black bars). These results indicate that the values computed by Bailer-Jones et al. (2018) cannot be used in our case; in general, distance values in Bailer-Jones et al. (2018) are useful in modelling the Galactic global structure by using many stars, but they are not well suited for analyzing individual halo stars when the parallax error is large.

Fig. 5. Distribution of Heliocentric distances. Results from Bailer-Jones et al. (2018) for sources in Muraveva et al. (2014), (top panel), and Pawlak et al. (2013), (bottom panel); the respective bin widths of the histograms are 0.38 kpc and 0.23 kpc (see the text for details). The vertical black bars signal the true values from Pietrzyński et al. (2013), (top panel), and Hilditch et al. (2005), (bottom panel). The cumulative relative frequency is plotted as a discontinuous curve.

In summary, we find strong statistical evidence for an affirmative answer to the question posed in the title of this section, at least for most of the sources in the sample selected in Sect. 2.3. This is particularly true for the 174 stars with median Galactocentric distance greater than 30 kpc, median Galactocentric Galactic velocity greater than 500 km s⁻¹, and well-behaved astrometric solutions.

2.6. Data processing pipeline validation

Our selection criteria are quite different from those used by Bromley et al. (2018), Hattori et al. (2018), or Marchetti et al. (2019) and we may have very few sources, if any, in common with them. Nonetheless and in order to assess the quality of our data processing techniques, we computed both Galactocentric distance and velocity for two stars discussed by Marchetti et al. (2019), Gaia DR2 4326973843264734208 and 4395399303719163904 (they are part of our 4 831 731 sources sample), in their Table 2 (see Table 1); our results for these stars are shown in Table 1 together with the results for Gaia DR2 4326973843264734208 obtained by Bromley et al. (2018) and they are statistically consistent with theirs. We also found that the angle between the position and velocity vectors of Gaia DR2 4326973843264734208 is 91°, which suggests an extragalactic origin if unbound as position and velocity directions are nearly perpendicular; this is also the interpretation favoured by Marchetti et al. (2019). As for Gaia DR2 4395399303719163904, we found an angle of 52°, which is compatible with an origin within the bulk of the Milky Way, but not at the centre, also matching the conclusion in Marchetti et al. (2019). Two other quality-control stars are shown in Table 1. Therefore, our data processing techniques are reliable enough to reproduce results previously obtained by other authors. In any case, the σ_π/π values of Gaia DR2 1383279090527227264, 4326973843264734208, 4395399303719163904, and 6456587609813249536 are 0.139, 0.155, 0.258, and 0.195, respectively; significantly below the typical values (>0.4) for the sources in our primary sample.

3. Full sample: results

Figure 6 shows the Galactocentric velocity components of the Gaia DR2 sample of 15 681 sources, our primary full sample. The components U, V, and W are positive in the directions of the Galactic centre, Galactic rotation, and the NGP. The sources in grey show good symmetry in the U − W plane but an obvious asymmetry in the U − V plane. Although the fraction of stars with V <  0 is 33.3%, their dispersion is wider than that of sources with V > 0 (meaning that there are more stars with extreme negative values of V), which we interpret as a signature of an unbound population as they move retrograde in V, this together with the large values of the velocity components. It is however worth to consider here that Schönrich et al. (2011) have shown in their Sect. 3.1 that, when the distance is uncertain, retrograde V components tend to be significantly larger than their prograde counterparts; in other words, the increased dispersion observed for V <  0 could be a side effect of the large uncertainties in d. On the other hand and in consistency with our initial hypothesis, 10 228 sources (∼65% of our primary full sample) have median Galactocentric speeds above 500 km s⁻¹, which is our threshold for hypervelocity candidacy (in other words, most distant stars are hypervelocity star candidates). If, instead of using the Monte Carlo methodology described above, the usual expressions are applied, 13 293 sources (∼85% of our primary full sample) have nominal Galactocentric speeds above 500 km s⁻¹.

Fig. 6. Galactocentric velocity components of the Gaia DR2 sample of 15 681 sources (in light grey). Sources with a ratio between the value of the parallax and its uncertainty >1 (1150) are plotted in black. The clustering centred at (U, V)∼(−50, −200) km s−1 corresponds to the Large Magellanic Cloud – (U, V, W)∼(−58, −221, 209) km s−1 – (mostly) and the Small Magellanic Cloud – (U, V, W)∼(26, −178, 175) km s−1 (values corrected for the Solar motion and the local standard of rest from those in Gaia Collaboration 2018b).

Unfortunately, the uncertainties associated with most sources are very large. Out of the primary full sample, we singled out 1150 sources (∼7%) for which the uncertainty in the value of the parallax is smaller than its nominal value in Gaia DR2. These sources are plotted in black in Fig. 6 and they exhibit a pattern of asymmetry consistent with the one found for the primary full sample. Figure 6, bottom panel, shows 745 sources with positive V-component (∼65%, these may have an origin in the Milky Way and be part of the halo) and 405 with negative values of V; among these, we observe sources that may be in the Magellanic Clouds – the clustering centred at (U, V)∼(−50, −200) km s⁻¹. The fact that the sources moving retrograde have higher net speeds suggests that they are indeed intergalactic stars although some substructure (but see the cautionary note above), other than that probably associated with the Magellanic Clouds, is also present in the form of kinematically coherent groups, perhaps debris from disrupted galaxies or actual dwarf galaxies in the process of falling towards the Milky Way galaxy or returning towards intergalactic space after experiencing a hyperbolic encounter with the Galaxy.

Figure 7 focusses on the best sample including the 174 stars, sources in black, with median Galactocentric distance greater than 30 kpc and median Galactocentric Galactic velocity greater than 500 km s⁻¹ that verify the criteria for reliable astrometry discussed in Sect. 2.3 and shown in Fig. 2. Out of 393 sources with reliable astrometric determination, 168 have V <  0; out of the 174 stars in the best sample (which are a subset of the 393), 112 have V <  0. Therefore, most stars in the best sample of 174 follow retrograde trajectories; the percentage of distant hypervelocity candidate stars among the astrometrically well-behaved subsample with V <  0 is 66.7%, which is the exact opposite of the trend found for the sample of 1150 sources for which the uncertainty in the value of the parallax is smaller than its nominal value in Gaia DR2 (Fig. 6, bottom panel).

Fig. 7. Same as Fig. 6 but focussing on the best sample including the 174 stars with median Galactocentric distance greater than 30 kpc and median Galactocentric Galactic velocity greater than 500 km s−1 that verify the criteria for reliable astrometry discussed in Sect. 2.3 and shown in Fig. 2. About 64% of stars in the best sample have negative median values of the V-component.

Figures 6 and 7, sources in black, suggest that there is indeed a population of fast and faraway stars that may have an extragalactic origin (mostly negative values of V for the best sample). This is consistent with the discussion in Cicone et al. (2016) and the results in Marchetti et al. (2019), and it appears to confirm the presence of a significant population of high-velocity stars of intergalactic provenance. The size of this unbound population must be very large because only the brightest (and therefore most massive) stars are present in the Gaia DR2 sample of very distant stars.

As for the spatial distribution of these distant hypervelocity star candidates, Figs. 8 and 9, top panels, show the location in the sky of our primary sample; the smaller samples, with lower uncertainties, are displayed in the panels at the bottom. The top panels show that the sources outline the Galactic disc, with most candidates observed projected towards the Galactic bulge (Galactic longitude, l <  50° and l > 310° in Fig. 9). There is however a non-negligible fraction of candidates faraway from the disc. Brown (2015) points out that the previously known hypervelocity stars exhibit statistically significant clustering towards the equatorial coordinates, α = 11^h30^m and δ = 3°. Such a concentration is not observed in our full sample and this could be a systematic effect or that the previously known hypervelocity stars belong to a separate population.

Fig. 8. Sources in equatorial coordinates. Distribution in equatorial coordinates of our primary sample of high-velocity and hypervelocity star candidates (top panel, only the sources with median Galactocentric Galactic velocity greater than 500 km s−1 are shown), those with a ratio between the value of the parallax and its uncertainty >1 (middle panel), and our best sample with most reliable astrometric solution (bottom panel). The colour map shows the values of the Galactocentric Galactic velocity in km s−1.

Fig. 9. Same as Fig. 8 but in Galactic coordinates.

What is observed in Fig. 9 is a scarcity of candidates towards the Galactic anticentre, l ∼ 180°, but this must be the result of a selection effect. Ejections from the region of the Galactic centre and bulge moving towards the anticentre will have large values of the radial velocity and comparatively negligible proper motions (many such sources may not have data in Gaia DR2). On the other hand, runaway exo-Galactic stars observed towards the Galactic anticentre may have virtually zero radial velocity if their paths do not cross that of the Galaxy; in other words, a value may be absent from Gaia DR2 and therefore they will not be part of our sample due to our selection criteria. The low number of candidates observed towards the Galactic anticentre could be a systematic effect. The presence of clustering is very obvious in the bottom panels of Figs. 8 and 9 and, without a doubt, this must be due to a systematic effect – some is expected towards the ecliptic poles, see the discussion in Lindegren et al. (2018), and the observed locations match approximately. The quality of the astrometric data acquired by the spacecraft for distant sources must be higher towards those directions for technical reasons. It cannot be due to the motion of the Milky Way galaxy – towards α = 10^h30^m and δ = −24°, with respect to the cosmic microwave background (Kogut et al. 1993) – as the Galactic apex (or antapex) is well separated from the concentrations of sources (see Fig. 8, top panel). The fact that most of the candidates appear projected towards the disc suggests that some may have an origin in it as discussed by Irrgang et al. (2018) and Marchetti et al. (2019). A number of sources appear projected towards the Large Magellanic Cloud (α = 5^h23^m, δ = −69 ° 45′, see Fig. 8, top panel) and they may have an origin in it, like HVS 3 (Irrgang et al. 2018; Erkal et al. 2019) and some B-type hypervelocity stars (Boubert et al. 2017); these sources are absent from the panels at the bottom due to the application of the various quality criteria discussed in Sect. 2.3.

4. Best samples: results

The number of high-velocity candidates, >500 km s⁻¹, in the sample including 1150 sources (for which the uncertainty in the value of the parallax is smaller than its nominal value) is 539 (∼47%); for the sample including 393 sources (those with the most reliable astrometric solutions) the number of high-velocity candidates is 174 (∼45%). In other words, for the sources with the lowest uncertainties the probability of being part of the high-velocity tail associated with the stellar halo or having a true hypervelocity nature is nearly the same.

Regarding the origin of the hypervelocity star candidates, Fig. 10 shows the distribution of Galactocentric distances and velocities as a function of the angle between the Galactocentric position and velocity vectors. Out of 539 sources, 293 have an angle <90°; out of 174, the equivalent number is 86. In addition, the mean value is 90°, the standard deviation is 6°, the median value is 90°, and the IQR is about 5° for the 539 sources (for the sample of 174, the respective results are 91°, 7°, 90°, and about 5°). These values suggest that the angular distribution is Gaussian. In any case, the most probable value of the angle between two randomly generated vectors is 90°. Nearly half of the best sources may have been ejected from the disc and the others may come from intergalactic space.

Fig. 10. Distances, velocities, and angles. Distribution of Galactocentric distances and velocities as a function of the angle between the Galactocentric position and velocity vectors for the full sample (in grey) of 15 681 sources and the one of 393 sources (in black) with the best astrometric solutions (left panels). The magnified version of the same data plot (right panels) includes the error bars for those sources with ratio between the value of the parallax and its uncertainty >2. Our lower cut-off limits, 500 km s−1 and 30 kpc, are indicated by a discontinuous line on each panel.

As for the Galactocentric speed, the average value is 786 km s⁻¹ with a standard deviation of 231 km s⁻¹, a median value of 737 km s⁻¹, and an IQR of 308 km s⁻¹ for the 539 sources (for the sample of 174, the respective results are 819 km s⁻¹, 246 km s⁻¹, 785 km s⁻¹, and 301 km s⁻¹). Although the overall uncertainties are large even for the sources with the best data (see the error bars in Fig. 10, right-hand panels), the statistical trends are clear and the 539 (174) sources can be considered hypervelocity star candidates although some contamination from the foreground stellar populations is expected. Figure 10 also shows that neither runaways from the Galactic centre nor stars moving towards it are present in the sample as these may be comparatively rare (but some high speed, under 500 km s⁻¹, candidates are present in Figs. 10 and 11). Figure 11 shows that a number of stars with V ∼ 250 km s⁻¹ have U ∼ −750 km s⁻¹ and θ ∼ 90°, which again hints at ejections from the Galactic disc.

Fig. 11. Velocity components and angles. Distribution of Galactocentric Galactic velocity components as a function of the angle between the Galactocentric position and velocity vectors for the full sample (in grey) of 15 681 sources and the one of 393 sources (in black) with the best astrometric solutions.

We found that 46 sources have uncertainties that are smaller than half the value of the distance, and 10 of them (see Table 2, top section) have values of the Galactocentric speed >500 km s⁻¹ (see Table 3, top section). Several of these sources appear to travel together, Gaia DR2 5303914835370054656, 5304019117178317312, 5317229302706969088, and 5317776481532378240, which is somewhat consistent with recent analyses by Hawkins & Wyse (2018) and Kamdar et al. (2019). The results for the 10 sources with the lowest uncertainties are still far from satisfactory in terms of errors. The uncertainty in distance (upper limit, meaning 84th percentile) is still close to its median value and the one in velocity is at least 75% its median value. The bottom section of Tables 2 and 3 shows three relevant examples with the highest values of the Galactocentric velocity and the lowest relative errors, which are at least 80%. All the stars in Table 3 have angles consistent with an origin in the disc (mostly) or outside the Milky Way galaxy. Table 2 shows that all the candidates have |b| < 23°, which is suggestive of an origin in the disc or, alternatively, a signal of possible contamination by disc stars with erroneous parallax determinations. The presence of a fraction of false positives among our candidates cannot be discarded (see an obvious case in Sect. 6), but any star leaving the disc tangent to its original bound trajectory may have an initially low value of its Galactic latitude (see also Figs. 10 and 11). For example, LAMOST-HVS4 (a hypervelocity star with a likely origin in the Galactic disc) has b = −17.92° (Li et al. 2018).

The vast majority of the candidates in Table 2 do not pass the astrometric quality filters discussed in Sect. 2.3. If we focus on the 174 sources, which have the best astrometric solutions, we obtain the best 10 candidates in Tables 4 and 5. Only one star, Gaia DR2 2130229806599268480, is also included in Tables 2 and 3. This shows that, for our data, having a well-behaved astrometric solution and σ_π/π <  0.5 are seldom compatible. The sample in Tables 4 and 5 includes one candidate with |b| >30°, Gaia DR2 1696697285206197248, which may have an origin in the Galactic bulge region (θ ∼ 68°). In addition, Gaia DR2 5348384273914768000 emerges as a clear candidate to have an exo-Galactic origin (θ ∼ 119°).

5. Statistical significance

Data from Gaia DR2 for sources beyond 30 kpc are affected by significant uncertainties, very large in the case of the parallax values (and therefore for the distances), although the values of the radial velocity and those of the components of the proper motion are sufficiently reliable in most cases. Given the large values of the uncertainties in the values of the Galactocentric distance and particularly in those of the velocity, even in the case of the most promising candidates (see Tables 2–5), our results must be interpreted with some caution. However, the fact that Fig. 6 shows a consistent asymmetry in the values of the V component of the Galactocentric velocity suggests that most of the sources in our samples could be unbound from the Milky Way galaxy. This is a statistically significant result that receives further support from the fact (see Sect. 2.5) that the net effect of the Eddington-Trumpler-Weaver bias makes the stars appear closer than they really are. If most of the sources in Tables 2–5 are more distant than estimated, their Galactocentric velocities are also probably higher; in other words, the median values of d_GC and v_GC in Tables 3 and 5 are mostly and probably lower limits and their true values could be close to the 84th percentile as shown in the tables.

It can be argued that our criteria to select hypervelocity candidates is perhaps too optimistic and that the statistical significance of our best candidates is somewhat marginal, however our approach is in line with the ones used in other studies; for example, one recent discovery by LAMOST has a Galactic rest-frame radial velocity of 408 km s⁻¹ at a Galactocentric radius of ∼30 kpc (Huang et al. 2017). Several candidates in Tables 3, 5, and 7 have comparable values even when the uncertainties are factored in. It is however true that LAMOST candidates are more reliable because their large velocities are virtually unaffected by the uncertainty in the stellar distance of their stars.

6. Discussion

Most of the distant sources with Galactocentric speed under 500 km s⁻¹ could be part of the regular halo of the Milky Way or be part of an outer spiral arm such as the one described by Dame & Thaddeus (2011). Star formation away from the Galactic disc is possible and it is a well-studied subject (for example, see de la Fuente Marcos & de la Fuente Marcos 2008a). Those in the halo may also have an external origin as the Milky Way cannibalizes nearby dwarf galaxies (Ibata et al. 1994). Du et al. (2018) find 24 high-velocity stars using Gaia DR2 data and spectroscopy from the LAMOST Data Release 5 and five of them appear to have an origin in the tidal debris of a disrupted dwarf galaxy, the rest seem to come from the Galactic disc. LAMOST-HVS1 has an intrinsic ejection velocity of 568 km s⁻¹ and it may come from a yet-to-be-discovered massive young star cluster located near the Norma spiral arm, although it has b = +35.41° (Hattori et al. 2019). LAMOST-HVS4 could be following a path with an origin in the disc (not the Galactic centre) with an outbound velocity as high as 590 ± 7 km s⁻¹ (Li et al. 2018). In general, sources with Galactocentric speed above 500 km s⁻¹ cannot be explained within the framework of synthetic models such as GALAXIA (Sharma et al. 2011), they must be the result of scattering or ejection events here in the Milky Way galaxy, within the Local Group, or elsewhere in the Universe.

Mechanisms capable of explaining the high-velocity stars coming from the disc include the dynamically driven double-degenerate double-detonation scenario (Shen et al. 2018), rapid mass-loss from a companion in a close binary including a supernova event (for example, see Tauris 2015) and interactions with intermediate-mass black holes within rich star clusters (Brown 2015; Boubert et al. 2018; Irrgang et al. 2018). The putative population of extragalactic origin may come from similar events outside our galaxy but also have more exotic sources such as disrupted galaxies, encounters with runaway supermassive black holes or intergalactic supernovae. Hypervelocity stars can be easily produced at the centres of galaxies via close encounters between binary stars and supermassive black holes; for example, HV 1 from the centre of the Milky Way (Brown et al. 2005) or HVS 3 from the centre of the Large Magellanic Cloud (Irrgang et al. 2018; Erkal et al. 2019), although a fly-by with a supermassive black hole may occur outside of a galaxy. Runaway supermassive black holes are expected to exist (Saslaw & De Young 1972; Saslaw et al. 1974; Mikkola & Valtonen 1990; Bonning et al. 2007; Gualandris & Merritt 2008) – but also intermediate-mass ones (Aarseth 2012) – and, in their way towards intergalactic space, they may scatter a significant number of stars. In addition, they may trigger star formation themselves on free-floating giant molecular clouds (de la Fuente Marcos & de la Fuente Marcos 2008b). Galaxy-less star formation may produce intergalactic supernovae (Gal-Yam et al. 2003), which in turn may produce runaway stars.

One of the strongest objections that one may make regarding the interpretation of our results is that the vast majority of candidates in Tables 2 and 4 have |b| < 23°. Tables 6 and 7 show candidates from the sample of 174 sources with well-behaved astrometric solutions that have |b| > 30°. One of them, Gaia DR2 1696697285206197248, was already included in Tables 4 and 5. Although in smaller numbers, candidates located well away from the Galactic disc are present in our samples as well (see also Fig. 9).

On the other hand, our sample is made of sources with estimated values of both the line-of-sight extinction and reddening in Gaia DR2; it is therefore possible to construct a colour-magnitude diagram (CMD) with the data. Figure 12, left-hand side panel, shows the resulting CMD: in grey, the full sample, in black, the 393 stars with the best astrometric solutions. This CMD has been obtained as the ones in Fig. 5 of Gaia Collaboration (2018c), Fig. 19 of Andrae et al. (2018), or Fig. 3 in de la Fuente Marcos & de la Fuente Marcos (2018). The median values were computed from the Monte Carlo sampling described above by assuming that the standard deviation in the estimated values of the line-of-sight extinction A_G and reddening E(G_BP − G_RP) from Gaia DR2 can be approximated by the difference between their respective 84th and 16th percentiles (which is a rather pessimistic assumption). The same approximation was used to estimate the error bars in Fig. 12, right-hand side panel. As a reference, five relevant PARSEC v1.2S + COLIBRI S_35 (Bressan et al. 2012; Marigo et al. 2017; Pastorelli et al. 2019)⁴ isochrones – of ages 1 Myr (in violet), 5 Myr (in purple), 10 Myr (in blue), 25 Myr (in cyan), and 50 Myr (in green) – are also plotted. The entries in Tables 2, 4, and 6 (see Fig. 12, right-hand side panel) appear to be consistent with being massive stars, younger than about 50 Myr. This young age together with their high velocities imply that they have been born in the Local Group, probably within 100 kpc from the centre of the Milky Way galaxy (unless they are blue stragglers). The fact that all the best candidates are consistently located within the CMD, in the region of the classical Cepheids, can be considered as a robust supporting argument in favour of the conclusions obtained in Sect. 2.5, where the effect of the Eddington-Trumpler-Weaver bias was discussed. Even sources out of the best sample exhibit conspicuous sequences that outline the various, temporarily stable, phases of stellar evolution for massive stars, including the loci of the luminous blue variables and yellow hypergiants with absolute magnitudes close to −10 mag.

Fig. 12. Colours, magnitudes, and theoretical isochrones. Colour-magnitude diagram (CMD) for sources in Fig. 6 (left panel), full sample of 15 681 sources (in grey) and the one of 393 sources (in black). Five isochrones of ages 1 Myr (in violet), 5 Myr (in purple), 10 Myr (in blue), 25 Myr (in cyan), and 50 Myr (in green), and solar metallicity are plotted as a reference (see the text for details). CMD for the candidates in Tables 2 (in red), 4 (in blue), and 6 (in brown) with error bars (right panel).

We have been unable to find any relevant data (such as observations by existing spectroscopic surveys) on most of the stars in Tables 2, 4, and 6, but future spectroscopic studies of some of the candidates found here may be able to confirm (or refute for false positives) an origin other than the Milky Way galaxy as discussed by Hawkins & Wyse (2018). The two exceptions are Gaia DR2 2079869356551980672 and 1360405705322797312. Gaia DR2 1360405705322797312 is a bona fide member of the metal poor and very old globular cluster NGC 6341 or M 92 as its parameters match well those of Cl* NGC 6341 BARN 73 (for example, see Pilachowski et al. 2000 ; Kamann et al. 2014); therefore, it is a false positive (likely induced by the crowded field) as the Galactocentric distance to NGC 6341 is 9.6 kpc, but the (wrong) value in Table 7 is 34.6 kpc.

In sharp contrast, the other previously studied object (meaning spectroscopically), Gaia DR2 2079869356551980672 or KIC 8828284, is given a spectral type of K5 II (bright giant) by Frasca et al. (2016) and K2 Ib (supergiant) by Gray et al. (2016) with V = 14.23 mag, which opens the door to a distance range from about 18 kpc (if bright giant with M_V = −2.0 mag) to perhaps 60 kpc (if supergiant with M_V = −4.7 mag). Our data processing pipeline (see Table 5) gives a value of the Galactocentric distance of 38.2 kpc for this object, which is consistent with the spectroscopic range, and a Galactocentric velocity of 869 km s⁻¹ for an angle of . It also gives and mag. In other words, if Gaia DR2 2079869356551980672 is a giant, it might still be bound to the Galaxy but if it is a supergiant, it is a robust example of a passing hypervelocity star with a velocity well above 1000 km s⁻¹ and an origin outside the Milky Way galaxy.

Regarding the other stars in Tables 2, 4, and 6, it may be argued that we have not discussed the reliability of these candidates in terms of the Re-normalised Unit Weight Error (RUWE) criterion, which is based on the chi squared of the astrometric fit and stellar colour (see technical note GAIA-C3-TN-LU-LL-124-01)⁵, and can be found from a number of Gaia DR2-related data servers⁶. Lindegren et al. (2018) argues that RUWE ≤ 1.4 is a reasonable criterion to identify sufficiently good solutions. The RUWE index is included in Tables 3, 5, and 7 and all the sources have RUWE < 1.4. Therefore, we confirm that our best candidates have well-behaved astrometric solutions.

7. Conclusions

In this paper, we have investigated the spatial distribution of distant (nominal Galactocentric distances greater than 30 kpc) hypervelocity star candidates present in Gaia DR2 data in an attempt to understand the origin of the putative anisotropic spatial distribution discussed by Brown et al. (2009) for example, and confirm or refute the presence of a significant population of high-velocity stars of intergalactic provenance as pointed out by Marchetti et al. (2019). No global parallax zero-point correction was performed. Our conclusions are:

• (i)

  The effect of the Eddington-Trumpler-Weaver bias is very important for samples such as the ones studied here. We estimate that over 80% of the sources discussed in our work may be located considerably further away than implied by the values of their parallaxes and, therefore, most of our velocity estimates are just lower limits.

• (ii)

  Most distant sources in Gaia DR2 have kinematics consistent with that of hypervelocity stars, but it cannot be discarded that some candidates may actually belong to the foreground stellar populations (like in the case of the false positive Gaia DR2 1360405705322797312, which is a member of the globular cluster NGC 6341).

• (iii)

  We have identified hundreds of distant hypervelocity star candidates that might be experiencing fly-bys to the Milky Way (their V-components are negative). The immediate neighbourhood of the Galaxy appears to host a significant population of exo-Galactic origin as pointed out by Marchetti et al. (2019) and such a population is consistent with findings by Cicone et al. (2016) for example, for a large sample of galaxies in the local Universe.

• (iv)

  We have been unable to identify a single hypervelocity star candidate with a probable origin in the region of the Galactic centre. The absence of such candidates is consistent with the results obtained by Irrgang et al. (2018) and Marchetti et al. (2019) for other samples of Gaia DR2 sources with shorter Galactocentric distances.

• (v)

  We have identified hundreds of distant high-velocity stars with a probable origin in the Galactic disc, confirming previous results by Boubert et al. (2018), Hattori et al. (2018, 2019), Irrgang et al. (2018), Li et al. (2018), and Marchetti et al. (2019).

• (vi)

  We confirm that the spatial distribution of hypervelocity star candidates is anisotropic, but the origin of such an anisotropy is probably the result of observational biases and selection effects.

• (vii)

  We identify a sample of candidate young, high-velocity stars with relatively low uncertainties that may have an intergalactic provenance. One of them, Gaia DR2 2079869356551980672 or KIC 8828284, emerges as a probable hypervelocity, supergiant star with a velocity well above 1000 km s⁻¹ and an origin outside the Milky Way galaxy.

• (viii)

  Most distant hypervelocity star candidates identified here have probable ages under 50 Myr. This young age together with the values of their velocities argue for birth places located within 100 kpc from the centre of the Milky Way.

As a closing argument, it is probably worth to consider that in this work we have put the stress on the consequences of the Eddington-Trumpler-Weaver bias, which leads to an underestimate of the heliocentric distance; however, we have to concede that the existence of a negative zero-point offset in the Gaia DR2 parallax has the opposite effect, leading to an overestimate of the heliocentric distance and, subsequently, the Galactocentric distance and velocities components. Which one, if any, is dominant in our case? We may not be able to find an answer to this difficult question in absolute terms (not even with the future Gaia DR3). It could be the case that the dominant effect (systematic overestimate due to Gaia parallax zero-point offset or underestimate due to Eddington-Trumpler-Weaver bias) may depend on the region of the sky under study or both effects may cancel each other out under certain circumstances or perhaps they follow some type of yet-to-be-found correlation. What appears to be clear is that the presence in our sample of reasonably robust hypervelocity star candidates such as Gaia DR2 2079869356551980672 (KIC 8828284) argues against concluding that Gaia DR2 is not of sufficient quality to fulfill our original science goals: performing a systematic exploration aimed at confirming or refuting that hypervelocity stars are ubiquitous in the local Universe and that, in the Milky Way, the known hypervelocity stars are anisotropically distributed.

Acknowledgments

We thank the referee for her/his constructive, actionable, and insightful reports that included very helpful suggestions regarding the presentation of this paper and the interpretation of our results, S. J. Aarseth for comments on black-hole-induced ejections and A. I. Gómez de Castro for comments on high-velocity stars and for providing access to computing facilities; RdlFM thanks L. Beitia-Antero for extensive discussions on Gaia DR2 data and F. J. Yáñez Gestoso for comments on some of the statistical issues discussed in this work. This research was partially supported by the Spanish “Ministerio de Economía y Competitividad” (MINECO) under grants ESP2015-68908-R and ESP2017-87813-R. In preparation of this paper, we made use of the NASA Astrophysics Data System, the ASTRO-PH e-print server, and the SIMBAD and VizieR databases operated at CDS, Strasbourg, France. 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.

References

Appendix A: Effects of a systematic parallax zero-point offset

As discussed in Sect. 2.2, there is admittedly robust evidence for the existence of a systematic negative parallax zero-point offset in Gaia DR2. This leads to values of the parallaxes that are too small, even after Monte Carlo sampling. Here, we used its most optimistic determination, −0.029 mas (Lindegren et al. 2018) that is based on quasars and was criticised by Schönrich et al. (2019) among others, to study its effects on our conclusions. In other words, we repeated the calculations described above but using d_c = 1/(π_c + 0.029 mas) instead of the usual d_c = 1/π_c. The immediate effect of this choice is that not a single one of the 174 high-velocity candidates (see Sect. 4) with the most reliable astrometric solutions survives as such after performing the parallax zero-point offset correction. However, 121 ± 2 of the 15 681 sources (nearly 0.8%) in our primary sample are still hypervelocity star candidates according to our criteria – at a Galactocentric distance >30 kpc, the Galactocentric velocity must be >500 km s⁻¹ – after performing the correction, but these are all sources with somewhat questionable astrometric solutions as they pass some quality-control criteria but not all of them (see Sect. 2.3). In sharp contrast and when no correction was applied, the fraction of hypervelocity star candidates is about 30% (4774 ± 8 out of 15681). In other words, 97% of the hypervelocity star candidates cannot be considered as such when the correction is performed. This fraction is comparable to the strength of the effect of the Eddington-Trumpler-Weaver bias discussed in Sect. 2.5, over 80% of the sources studied in our work may be located further away than implied by the values of their parallaxes.

If instead of using −0.029 mas, we consider the most pessimistic determination, −0.082 ± 0.033 mas (Stassun & Torres 2018), then the number of high-velocity candidates goes down to zero. The same result is obtained if we apply the zero-point offset computed by Schönrich et al. (2019), −0.054 ± 0.006 mas. This outcome was expected because the 87 733 672 sources with estimated values of the line-of-sight extinction and reddening in Gaia DR2 have strictly positive values of the parallax and these values of the zero-point offset (in absolute terms) are above the cutoff value 0.033 mas (or 30 kpc). This fact also explains why a small fraction of high-velocity candidates manage to survive the correction process when the most optimistic determination, −0.029 mas, is used.

Therefore, if one believes that not having carried out the zero-point offset correction is a serious weakness in our analysis, then one must conclude that Gaia DR2 may not be not adequate to fulfill our original science goals, perhaps the future Gaia DR3 will be better suited for the task, but see our closing argument in Sect. 7.

All Tables

All Figures

	Fig. 1. Analytical versus Monte Carlo distances for 100 representative synthetic data with σπ/π ∈ [0.2, 4.2] and π ∈ [0.0125, 0.0333] mas (uniformly distributed). These are shown with (left) and without (right) error bars, computed as described in the text. The discontinuous line shows the diagonal (d = dc) for reference.
In the text

Fig. 2. Quality summary of the primary sample (15 681 sources). Position and relative parallax error as a function of the parameter astrometric_excess_noise_sig or significance of excess noise (left panels). Lindegren et al. (2018) indicate that astrometrically well-behaved sources have a significance of excess noise ≤2 (region on the left of the red line). The green line in the left-hand, bottom panel signals σπ/π = 0.2. Dependency between quality flags (right panels). Lindegren et al. (2018) state that sources with astrometric_gof_al or goodness of fit <3 (below the green line of the top panel), mean_varpi_factor_al or mean parallax factor AL in the interval [−0.23, 0.32] (between the green lines of the middle panel), and visibility_periods_used or visibility periods >8 (above the green line of the bottom panel) should be used. Not shown here, but also applied, is rv_nb_transits or number of observations used to compute the value of the radial velocity >5.

In the text

Fig. 3. Extinction and reddening, a comparison. Left panels: extinction values AV from Schlegel et al. (1998), (top panels), and Schlafly & Finkbeiner (2011), (bottom panels), versus the line-of-sight extinction AG from Gaia DR2 for the 15 681 sources in our primary sample (in grey) and the 393 stars with the best astrometric solutions (in black). Right panels: colour excesses E(B − V) from Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) versus E(GBP − GRP) from Gaia DR2 for the same sources. The reference line in red is the diagonal. See the text for additional details.

In the text

Fig. 4. Eddington-Trumpler-Weaver bias. Evaluation of this bias for data from Muraveva et al. (2014), (top panel), and Pawlak et al. (2013), (bottom panel). The y-axis shows the difference between the value of the parallax from Gaia DR2 and the one derived from the assumed distance – 0.0208 ± 0.0005 mas for the LMC (top panel, Pietrzyński et al. 2013) and 0.0165 ± 0.0010 mas for the SMC (bottom panel, Hilditch et al. 2005)— divided by the quoted uncertainty (0.0005 mas, top panel, and 0.0010 mas, bottom panel). The x-axis shows the corresponding value of the relative error in parallax from Gaia DR2. Only sources with strictly positive values of the parallax are plotted; a sample of the full catalogues is shown as the full ranges in x and y are wider. The condition σπ/π = 0.2 is plotted as a discontinuous line. Left panels: results without considering any global zero-point correction; right panels: include the correction (π + 0.029 mas) discussed in the text.

In the text

Fig. 5. Distribution of Heliocentric distances. Results from Bailer-Jones et al. (2018) for sources in Muraveva et al. (2014), (top panel), and Pawlak et al. (2013), (bottom panel); the respective bin widths of the histograms are 0.38 kpc and 0.23 kpc (see the text for details). The vertical black bars signal the true values from Pietrzyński et al. (2013), (top panel), and Hilditch et al. (2005), (bottom panel). The cumulative relative frequency is plotted as a discontinuous curve.

In the text

Fig. 6. Galactocentric velocity components of the Gaia DR2 sample of 15 681 sources (in light grey). Sources with a ratio between the value of the parallax and its uncertainty >1 (1150) are plotted in black. The clustering centred at (U, V)∼(−50, −200) km s−1 corresponds to the Large Magellanic Cloud – (U, V, W)∼(−58, −221, 209) km s−1 – (mostly) and the Small Magellanic Cloud – (U, V, W)∼(26, −178, 175) km s−1 (values corrected for the Solar motion and the local standard of rest from those in Gaia Collaboration 2018b).

In the text

	Fig. 7. Same as Fig. 6 but focussing on the best sample including the 174 stars with median Galactocentric distance greater than 30 kpc and median Galactocentric Galactic velocity greater than 500 km s−1 that verify the criteria for reliable astrometry discussed in Sect. 2.3 and shown in Fig. 2. About 64% of stars in the best sample have negative median values of the V-component.
In the text

Fig. 8. Sources in equatorial coordinates. Distribution in equatorial coordinates of our primary sample of high-velocity and hypervelocity star candidates (top panel, only the sources with median Galactocentric Galactic velocity greater than 500 km s−1 are shown), those with a ratio between the value of the parallax and its uncertainty >1 (middle panel), and our best sample with most reliable astrometric solution (bottom panel). The colour map shows the values of the Galactocentric Galactic velocity in km s−1.

In the text

	Fig. 9. Same as Fig. 8 but in Galactic coordinates.
In the text

Fig. 10. Distances, velocities, and angles. Distribution of Galactocentric distances and velocities as a function of the angle between the Galactocentric position and velocity vectors for the full sample (in grey) of 15 681 sources and the one of 393 sources (in black) with the best astrometric solutions (left panels). The magnified version of the same data plot (right panels) includes the error bars for those sources with ratio between the value of the parallax and its uncertainty >2. Our lower cut-off limits, 500 km s−1 and 30 kpc, are indicated by a discontinuous line on each panel.

In the text

	Fig. 11. Velocity components and angles. Distribution of Galactocentric Galactic velocity components as a function of the angle between the Galactocentric position and velocity vectors for the full sample (in grey) of 15 681 sources and the one of 393 sources (in black) with the best astrometric solutions.
In the text

Fig. 12. Colours, magnitudes, and theoretical isochrones. Colour-magnitude diagram (CMD) for sources in Fig. 6 (left panel), full sample of 15 681 sources (in grey) and the one of 393 sources (in black). Five isochrones of ages 1 Myr (in violet), 5 Myr (in purple), 10 Myr (in blue), 25 Myr (in cyan), and 50 Myr (in green), and solar metallicity are plotted as a reference (see the text for details). CMD for the candidates in Tables 2 (in red), 4 (in blue), and 6 (in brown) with error bars (right panel).

In the text