Gaia DR2 in 6D: Searching for the fastest stars in the Galaxy

T. Marchetti, E. M. Rossi, A. G. A. Brown

Introduction

Stars with extremely high velocities have been long studied to probe our Galaxy. The interest in the high velocity tail of the total velocity distribution of stars in our Milky Way is twofold. First, it flags the presence of extreme dynamical and astrophysical processes, especially when the velocity of a star is so high that it approaches (or even exceeds) the escape speed from the Galaxy at its position. Secondly, high velocity stars, spanning a large range of distances, can be used as dynamical tracers of integral properties of the Galaxy. The stellar high velocity distribution has for example been used to trace the local Galactic escape speed and the mass of the Milky Way (Smith et al. 2007; Gnedin et al. 2010; Piffl et al. 2014, e.g.). To put the concept of high velocity in context, the value of the escape speed is found to be $\sim 530$ km s-1 at the Sun position, it increases up to $\sim 600$ km s-1 in the central regions of the Galaxy, and then falls down to $\lesssim 400$ km s-1 at Galactocentric distances $\sim 50$ kpc (Williams et al. 2017).

A first class of objects that can be found in the high tail of the total velocity distribution is fast halo stars. Their measured dispersion velocity is around 150 km s-1 (Smith et al. 2009; Evans et al. 2016), therefore 3- $\sigma$ outliers can exceed $450$ km s-1 , while remaining bound. Halo stars could also reach unbound velocities, when they are part of the debris of tidally disrupted satellite galaxies, like the Sagittarius Dwarf galaxy, that has not yet virialized (Abadi et al. 2009, e.g.). Velocities outliers in the bulge and disk velocity distribution may also exist and become apparent in a large data set.

"Runaway stars" (RSs) form an another class of high velocity stars. They were originally introduced as O and B type stars ejected from the Galactic disk with velocities higher than $40$ km s-1 (Blaauw 1961). Theoretically, there are two main formation channels: i) dynamical encounters between stars in dense stellar systems such as young star clusters (Poveda et al. 1967; Leonard & Duncan 1990; Gvaramadze et al. 2009, e.g.), and ii) supernova explosions in stellar binary systems (Blaauw 1961; Portegies Zwart 2000, e.g.). Both mechanisms have been shown to occur in our Galaxy (Hoogerwerf et al. 2001). Typical velocities attained by the two formation channels are of the order of a few tens of km s-1 , and even if several hundreds of km s-1 can be attained for the most extreme systems (Portegies Zwart 2000; Przybilla et al. 2008; Gvaramadze et al. 2009; Gvaramadze & Gualandris 2011; Silva & Napiwotzki 2011), simulations indicate that the majority of runaway stars from dynamical encounters have ejection velocities $\lesssim 200$ km s-1 (Perets & Šubr 2012). Recent results show that it is possible to achieve ejection velocities up to $\sim 1300$ km s-1 for low-mass G/K type stars in very compact binaries (Tauris 2015). Nevertheless, the rate of production of unbound RSs, referred to as hyper runaway stars (HRSs), is estimated to be as low as $8\cdot 10^{-7}$ yr-1 (Perets & Šubr 2012; Brown 2015).

As a class, the fastest stars in our Galaxy are expected to be hypervelocity stars (HVSs). These were first theoretically predicted by Hills 1988 as the result of a three-body interaction between a binary star and the massive black hole in the Galactic Centre (GC), Sagittarius A∗. Following this close encounter, a star can be ejected with a velocity $\sim 1000$ km s-1 , sufficiently high to escape from the gravitational field of the Milky Way (Kenyon et al. 2008; Brown 2015). The first HVS candidate was discovered by Brown et al. 2005: a B-type star with a velocity more than twice the Galactic escape speed at its position. Currently about $\sim 20$ unbound HVSs with velocities $\sim 300$ - $700$ km s-1 have been discovered by targeting young stars in the outer halo of the Milky Way (Brown et al. 2014). In addition, tens of mostly bound candidates have been found at smaller distances but uncertainties prevent the precise identification of the GC as their ejection location (Hawkins et al. 2015; Vickers et al. 2015; Zhang et al. 2016; Marchetti et al. 2017; Ziegerer et al. 2017, e.g.). HVSs are predicted to be ejected from the GC with an uncertain rate around $10^{-4}$ yr-1 (Yu & Tremaine 2003; Zhang et al. 2013), two orders of magnitude larger than the rate of ejection of runaway stars with comparable velocities from the stellar disk (Brown 2015). Because of their extremely high velocities, HVS trajectories span a large range of distances, from the GC to the outer halo. Thus HVSs have been proposed as tools to study the matter distribution in our Galaxy (Gnedin et al. 2005; Sesana et al. 2007; Kenyon et al. 2014; Rossi et al. 2017; Fragione & Loeb 2017; Contigiani et al. 2018, e.g.) and the GC environment (Zhang et al. 2013; Madigan et al. 2014, e.g.), but a larger and less observationally biased sample is needed in order to break degeneracies between the GC binary content and the Galactic potential parameters (Rossi et al. 2017). Using the fact that their angular momentum should be very close to zero, HVSs have also been proposed as tools to constrain the Solar position and velocity (Hattori et al. 2018a). Other possible alternative mechanisms leading to the acceleration of HVSs are the encounter between a single star and a massive black hole binary in the GC (Yu & Tremaine 2003; Sesana et al. 2006; Sesana et al. 2008, e.g.), the interaction between a globular cluster with a single or a binary massive black hole in the GC (Capuzzo-Dolcetta & Fragione 2015; Fragione & Capuzzo-Dolcetta 2016), and the tidal interaction of a dwarf galaxy near the center of the Galaxy (Abadi et al. 2009). Another possible ejection origin for HVSs and high velocity stars in our Galaxy is the Large Magellanic Cloud (Boubert & Evans 2016; Boubert et al. 2017; Erkal et al. 2018, LMC, ), orbiting the Milky Way with a velocity $\sim 380$ km s-1 (van der Marel & Kallivayalil 2014).

In addition to the unbound population of HVSs, all the ejection mechanisms mentioned above predict also a population of bound HVSs (BHVSs): stars sharing the same formation scenario as HVSs, but with an ejection velocity which is not sufficiently high to escape from the whole Milky Way (Bromley et al. 2006, e.g.). Most of the deceleration occurs in the inner few kpc due to the bulge potential (Kenyon et al. 2008), and the minimum velocity necessary at ejection to be unbound is of the order of $\sim 800$ km s-1 (Brown 2015; Rossi et al. 2017, a precise value depends on the choice of the Galactic potential,). If we consider the Hills mechanism , this population of bound stars is expected to be dominant over the sample of HVSs (Rossi et al. 2014; Marchetti et al. 2018).

At the moment, the fastest star discovered in our Galaxy is US 708, traveling away from the Milky Way with a total velocity $\sim 1200$ km s-1 (Hirsch et al. 2005). Its orbit is not consistent with coming from the GC (Brown et al. 2015), and the most likely mechanism responsible for its acceleration is the explosion of a thermonuclear supernova in an ultra-compact binary in the Galactic disk (Geier et al. 2015).

Using Gaia DR2 data, Boubert et al. 2018 show that almost all the previously discovered late-type HVS candidates are most likely bound to the Galaxy, and their total velocity was previously overestimated because of inaccurate parallaxes and/or proper motions. Only one late-type star, LAMOST J115209.12+120258.0 (Li et al. 2015), is most likely unbound, but the Hills mechanisms is ruled out as a possible explanation of its extremely high velocity. The majority of B-type HVSs from (Brown et al. 2014; Brown et al. 2015) are still found to be consistent with coming from the GC when using Gaia DR2 proper motions (Erkal et al. 2018).

In this paper we search for the fastest stars in the Milky Way, within the sample of $\sim 7$ million stars with a six-dimensional phase space measurement in Gaia DR2. Since the origin of high velocity stars in our Galaxy is still a puzzling open question, we simply construct the total velocity distribution in the Galactic rest-frame in order to identify and characterize the high velocity tail. In doing so, we do not bias our search towards any specific class of high velocity stars.

This manuscript is organized as follows. In Section 2, we explain how we determine distances and total velocities in the Galactic rest frame for the whole sample of stars. We presents results in terms of stellar total velocity in Section 3. In Section 4, we focus on the high velocity stars in the sample, and then in Section 5 we concentrate on the stars with a probability greater than $80\%$ of being unbound from the Galaxy, discussing individually the most interesting candidates. Finally, we conclude and discuss our results and findings in Section 6.

Distance and Total Velocity Determination

The Gaia catalogue provides parallaxes, and thus a conversion to a distance is required to convert the apparent motion of an object on the celestial sphere to a physical motion in space, that is needed to determine the total velocity of a star. Bailer-Jones 2015 discusses in details how this operation is not trivial when the relative error in parallax, $f\equiv\sigma_{\varpi}/\varpi$ , is either above $20\%$ or it is negative. We choose to separate the discussion on how we determine distances and total velocities of stars with $0<f\leq 0.1$ (the "low-f sample") and of those with either $f>0.1$ or $f<0$ (the "high-f sample"). There are $7183262$ stars with both radial velocity and the astrometric parameters (parallax and proper motions) in Gaia DR2, therefore in the following we will focus on this subsample of stars.

$5393495$ out of $7183262$ stars ( $\sim 75\%$ ) with radial velocity measurement in Gaia DR2 have a relative error in parallax $0<f\leq 0.1$ . For this majority of stars we can get an accurate determination of their distance just by inverting the parallax: $d=1/\varpi$ (Bailer-Jones 2015). We then model the proper motions and parallax distribution as a multivariate Gaussian with mean vector:

2 The "high-f Sample"

A more careful analysis is required for $1789767$ stars ( $\sim 25\%$ ) with either $f>0.1$ or with a negative measured parallax. For these stars, we follow the approach outlined in Bailer-Jones 2015; Astraatmadja & Bailer-Jones 2016a; Astraatmadja & Bailer-Jones 2016b; Luri et al. 2018; Bailer-Jones et al. 2018. We use a full Bayesian analysis to determine the posterior probability $P(d|\varpi,\sigma_{\varpi})$ of observing a star at a distance $d$ , given the measured parallax $\varpi$ and its Gaussian uncertainty $\sigma_{\varpi}$ . The authors show how the choice of the prior probability on distance $P(d)$ can seriously affect the shape of the posterior distribution, and therefore lead to significantly different values for the total velocity of a star. We decide to adopt an exponentially decreasing prior:

which has been shown to perform best for stars further out than $\sim 2$ kpc (Astraatmadja & Bailer-Jones 2016b), that is the expected distance of stars with a large relative error on parallax (see Appendix A). The value of the scale length parameter $L$ is fixed to $2600$ pc, and we refer the reader to the discussion in Appendix A for the reasons behind our choice of this particular value. By means of Bayes’ theorem we can then express the posterior distribution on distances as:

where the likelihood probability $P(\varpi|d,\sigma_{\varpi})$ is a Gaussian distribution centered on $1/d$ :

In our case, we decide to fully include the covariance matrix between the astrometric properties, following the approach introduced in Marchetti et al. 2017. In this case, for each star the likelihood probability is a three dimensional multivariate Gaussian distribution with mean vector:

and covariance matrix given by equation (2). The prior distribution on distance is given by equation (3), and we assume uniform priors on proper motions. We then draw proper motions and distances from the resulting posterior distribution using the affine invariant ensemble Markov chain Monte Carlo (MCMC) sampler emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013). We run each chain using $32$ walkers and $100$ steps, for a total of $3200$ random samples drawn from the posterior distribution. We initialize the walkers to random positions around the mean value of the proper motions and of the inverse of the mode of the posterior distribution in distance, equation (4), to achieve a fast convergence of the chain. We run $500$ burn-in steps to let the walkers explore the parameter space, and then we use the final positions as initial conditions for the proper MC chain. We then directly use this MC sampling to derive a distribution for the total velocity in the Galactic rest frame of each star, assuming the same parameters for the Sun presented in Section 2.1. We check that the mean acceptance fraction (i.e. the fraction of steps accepted for each walker) is between $0.25$ and $0.5$ as a test for the convergence of each MC chain (Foreman-Mackey et al. 2013).

The Total Velocity Distribution of Stars in Gaia DR2

Using the approach discussed in Section 2, we publish a catalogue with distances and velocities in the Galactocentric frame for all the $7183262$ stars analyzed in this paper. This is publicly available at http://home.strw.leidenuniv.nl/˜marchetti/research.html. A full description of the catalogue content can be found in Appendix B.

Figure 3 shows the Toomre diagram for all the $\sim 7$ million stars, a plot that is useful to distinguish stellar populations based on their kinematics. On the $x$ -axis we plot the component $V$ of the Galactocentric Cartesian velocity, and on the $y$ -axis the component orthogonal to it, $\sqrt{U^{2}+W^{2}}$ . Not surprisingly, most of the stars behave kinematically as disk stars on rotation-supported orbits, with $V$ values around the Sun’s orbital velocity (Gaia Collaboration et al. 2018b, see). A sub-dominant, more diffuse, population of stars with halo-like kinematics is also present, centered around $V=0$ and with a larger spread in total velocity.

High Velocity Stars in Gaia DR2

$-0.23\leq$ mean_varpi_factor_al $\leq 0.32$ ;

In order to get hints on the ejection location of our sample of high velocity stars, we perform numerical orbit integration of their trajectories back in time using the python package Gala (Price-Whelan 2017). For each star we use $1000$ random samples from the proper motions, distance, and radial velocity MC sampling discussed in Section 2. We integrate each orbit back in time for a total time of $1$ Gyr, with a fixed time-step of $0.1$ Myr, using the gala potential MilkyWayPotential. This is a four components Galactic potential model consisting of a Hernquist bulge and nucleus (Hernquist 1990):

where $i=b,n$ for the bulge and the nucleus, respectively, a Miyamoto-Nagai disk (Miyamoto & Nagai 1975):

and a Navarro-Frenk-White halo (Navarro et al. 1996):

The parameters are chosen to fit the enclosed mass profile of the Milky Way (Bovy 2015), and are summarized in Table 1. We then derive the pericenter distance and, for bound MC realizations, the apocenter distance and the eccentricity of the orbit. We also record the energy and the angular momentum of each MC orbit. We check for energy conservation as a test of the accuracy of the numerical integration.

Unbound Stars: Hypervelocity and Hyper Runaway Star Candidates

We find that all of these $7$ stars have orbits that, when integrated back in time, are not consistent with coming from the GC. Therefore, according to our classification criterion, there are no stars classified as HVS candidates. The absence of HVS candidates in the subset of Gaia DR2 with radial velocities was anticipated by predictions by Marchetti et al. 2018, analyzing the Hills mock catalogue of HVSs. This is due to the fact that the expected number density of HVSs generated via the Hills’ mechanism is expected to increase linearly with increasing galactocentric distance (Brown 2015), and the majority of HVSs in the Milky Way are too faint to have a radial velocity measurement from Gaia DR2. We cannot exclude the presence of bound HVSs in the subset of $\sim 7$ million stars considered in this work, but their identification is not trivial because of their complex orbits and lower velocities. About 20 BHVSs are expected to have radial velocities from Gaia DR2 (Marchetti et al. 2018), but their identification is beyond the scope of this manuscript.

2 Extragalactic Stars

Conclusions

None of these 20 stars is consistent with coming from the inner $1$ kpc, so there are no HVS candidates. This is consistent with estimates presented in Marchetti et al. 2018.

$13$ out of the $20$ unbound candidates have probabilities $<50\%$ to originate from the stellar disk of the Galaxy. This surprising and unexpected population of stars could be either produced as RSs/HRSs/HVSs from the LMC, thanks to its high orbital velocity around the Milky Way, or could be members of dwarf galaxies tidally disrupted by the gravitational interaction with the Galaxy. Further analyses are required in order to identify their origin.

This paper is just a first proof of the exciting discoveries that can be made mining the Gaia DR2 catalogue. We only limited our search to the $\sim 7$ million stars with a full phase space information, a small catalogue compared to the full $1.3$ billion sources with proper motions and parallaxes. Synergies with existing and upcoming ground-based spectroscopic surveys will be essential to obtain radial velocities and stellar spectra for subsets of these stars (Dalton 2016; de Jong et al. 2016; Kunder et al. 2017; Martell et al. 2017, e.g.). For what concerns HVSs, Marchetti et al. 2018 shows how the majority of HVSs expected to be found in the Gaia catalogue are actually fainter than the limiting magnitude for radial velocities in DR2. We therefore did not expect to discover the bulk of the HVS population with the method outlined in this paper, but other data mining techniques need to be implemented in order to identify them among the dominant background of bound, low velocity stars (Marchetti et al. 2017, see for example). We also show how particular attention needs to be paid to efficiently filter out contaminants and instrumental artifacts, which might mimic high velocity stars at a first inspection.

Acknowledgements

We thank the anonymous referee for his/her comments, which greatly improved the quality of this manuscript. We also thank E. Zari and the Gaia group meeting at Leiden Observatory for useful comments, suggestions and discussions during the preparation of this paper. TM and EMR acknowledge support from NWO TOP grant Module 2, project number 614.001.401. This project was developed in part at the 2017 Heidelberg Gaia Sprint, hosted by the Max-Planck-Institut für Astronomie, Heidelberg. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research made use of Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2013). All figures in the paper were produced using matplotlib (Hunter 2007) and Topcat (Taylor 2005). This work would not have been possible without the countless hours put in by members of the open-source community all around the world.

References

Appendix A Choice of the Prior Probability on Distances

We find $352010$ of the $7$ million stars to have $f=\sigma_{\varpi}/\varpi>0.1$ . We can see that this value is about 5 times smaller than the one found in Gaia DR2 (see Section 2.2). All these stars are found at distance larger than $\sim 4.5$ kpc from the Sun, and therefore we choose to adopt the exponentially decreasing prior to derive their distances (Astraatmadja & Bailer-Jones 2016b), see equation (3). The mode of the posterior distribution in equation (4) can be determined by numerically finding the roots of the implicit equation (Bailer-Jones 2015):

The reason why we choose not to use distances from Bailer-Jones et al. 2018 is that the authors fit the values of the scale length $L$ to a full three-dimensional model of the Galaxy Note that Bailer-Jones et al. 2018 adopt a scale length that varies smoothly with Galactic longitude and latitude.. Their values are therefore driven by nearby, bright disk stars, with $f\ll 1$ . Such an approach would underestimate distances (and therefore total velocities) to faint distant stars, the ones we are more interested in.

Appendix B Content of the Distance and Velocity Catalogue

Table 3 provides an explanation of the content of the catalogue containing distances and velocities for the $7183262$ stars with a radial velocity measurement in Gaia DR2. The catalogue is publicly available at http://home.strw.leidenuniv.nl/˜marchetti/research.html.