New candidate hypervelocity red clump stars in the inner Galactic bulge Open Access

ABSTRACT

We search for high-velocity stars in the inner region of the Galactic bulge using a selected sample of red clump stars. Some of those stars might be considered hypervelocity stars (HVSs). Even though the HVSs ejection relies on an interaction with the supermassive black hole (SMBH) at the centre of the Galaxy, there are no confirmed detections of HVSs in the inner region of our Galaxy. With the detection of HVSs, ejection mechanism models can be constrained by exploring the stellar dynamics in the Galactic centre through a recent stellar interaction with the SMBH. Based on a previously developed methodology by our group, we searched with a sample of preliminary data from version 2 of the Vista Variables in the Via Lactea (VVV) Infrared Astrometric Catalogue (VIRAC2) and Gaia DR3 data, including accurate optical and near-infrared proper motions. This search resulted in a sample of 46 stars with transverse velocities larger than the local escape velocity within the Galactic bulge, of which four are prime candidate HVSs with high-proper motions consistent with being ejections from the Galactic centre. Adding to that, we studied a sample of reddened stars without a Gaia DR3 counterpart and found 481 stars with transverse velocities larger than the local escape velocity, from which 65 stars have proper motions pointing out of the Galactic centre and are candidate HVSs. In total, we found 69 candidate HVSs pointing away from the Galactic centre with transverse velocities larger than the local escape velocity.

1 INTRODUCTION

The Galactic bulge is the region 2.5 kpc around the Galactic centre (GC), and within the Λ cold dark matter paradigm, it was the first Galactic structure that formed through the hierarchical merging of less-massive structures (Steinmetz & Muller 1995; Springel & Hernquist 2005; Barbuy, Chiappini & Gerhard 2018). The environment in such region is extreme, with high stellar crowding and in its centre resides a supermassive black hole (SMBH), Sagittarius A* (Sgr A*), with a mass of |$4.3\times 10^6\, \mathrm{{\rm M}_{\odot }}$| (Rodriguez 1978; Ghez et al. 2008; Genzel, Eisenhauer & Gillessen 2010; GRAVITY Collaboration 2022). The high stellar density favours interactions amongst stars and between stars and Sgr A*.

An interaction between a triple system can result in one of the stars acquiring a velocity larger than the local standard of rest. The identification of a star after the interaction gives insight into the dynamics surrounding the interaction itself. For example Blaauw (1961) showed that runaway stars with peculiar velocities larger than 30 |$\mathrm{km\, s^{-1}}$| can form after a binary stellar system disruption in the supernova explosion of one of its components. Poveda, Ruiz & Allen (1967) investigated dynamical stellar interactions during the collapse of small clusters of massive stars, noting that runaway stars can result from such interactions acquiring velocities that exceed 35 |$\mathrm{km\, s^{-1}}$| and in some cases reaching nearly 200 |$\mathrm{km\, s^{-1}}$|⁠.

Hypervelocity stars (HVSs) are the fastest stars in the Galaxy. They were initially defined as stars ejected after a three-body interaction of a binary system with Sgr A*; one star is ejected as an HVS, while the other remains attached to the SMBH as an S-Star¹ (Hills 1988). This is known as the Hills mechanism, the most accepted and successful HVS ejection scenario, but not the only one. Alternative mechanisms that explain how an HVS acquires its large velocity include the interaction between a binary massive black hole (bMBH) and a single star (Yu & Tremaine 2003); or between a globular cluster and an SMBH (Brown 2015; Capuzzo-Dolcetta & Fragione 2015; Fragione, Capuzzo-Dolcetta & Kroupa 2017; Irrgang et al. 2019; Neunteufel 2020). In a tidal disruption event, the ejection velocity of a binary component is proportional to the binary components separation and the mass of the BH: v_ej ∝ (a)^−1/2(M/M_⊙)^1/6. Given the mass of Sgr A*, |$4\times 10^6\, \mathrm{{\rm M}_{\odot }}$|⁠, HVSs can be ejected at up to |$\sim 4000$||$\mathrm{km\, s^{-1}}$| (Hills 1988; Sari, Kobayashi & Rossi 2010; Rossi et al. 2017, 2021), which is larger than the escape velocity. The local escape velocity is ∼830 |$\mathrm{km\, s^{-1}}$| in the GC, ∼320 |$\mathrm{km\, s^{-1}}$| in the halo (50 kpc), and ∼530 |$\mathrm{km\, s^{-1}}$| at the Sun location (Rossi et al. 2017; Deason et al. 2019). Besides explaining the extreme velocities of HVSs, the Hills mechanism can explain the presence of young stars close to the GC, where an in situ formation is unlikely because a molecular cloud would not survive such an extreme environment. In this scenario, the young stars close to the GC could be the remnants of a binary system disruption by Sgr A* (Hills 1988; Generozov 2021).

Since the first discovery by Brown et al. (2005), there are around 20 confirmed HVSs and over 500 candidates (Boubert et al. 2018). The high stellar crowding as well as large and patchy extinction in the central regions of the Milky Way hampered the detection of HVSs close to their origin. They were identified in a more favourable environment, the Milky Way halo, with velocities ranging from 300 up to 1700 |$\mathrm{km\, s^{-1}}$|⁠, exceeding the local escape speed (e.g. Kollmeier et al. 2009; Kenyon et al. 2014; Palladino et al. 2014; Brown et al. 2018). Almost all of them are B-type stars with masses between 2.5 and 4 M_⊙, except for the fastest star amongst them: HVS-S5 (Koposov et al. 2020), an A-type star with a velocity of 1700 |$\mathrm{km\, s^{-1}}$|⁠. This star is also the only HVS whose travel direction and orbit points to the origin from the GC, favouring the Hills mechanism as responsible for its ejection. For other cases, the Hills mechanism can be ruled out, as the orbits of the stars suggest an origin in the disc, in the Magellanic Clouds (e.g. Przybilla et al. 2008; Lennon et al. 2017; Boubert, Erkal & Gualandris 2020; Evans et al. 2021), or the Sagittarius dwarf spheroidal galaxy (Huang et al. 2021; Li et al. 2022).

In a recent study, Li et al. (2023) combined multiple spectroscopic large surveys and found that the ejection from dwarf galaxies and globular clusters is more prominent for the metal-poor late-type halo HVS production. Based on Gaia data, and follow-up spectroscopic observations, El-Badry et al. (2023) identified four hypervelocity white dwarf (WD) stars, travelling at space velocities larger than 1300 |$\mathrm{km\, s^{-1}}$| after the double detonation – helium in the surface, followed by carbon in the core – of the more massive WD in a binary system of WDs. The objects have crossed the Galactic plane on multiple occasions, getting accelerated or slowed down, but none of them have trajectories that trace back to the GC. The analysis of positions and velocity distributions of the fastest stars in proximity to the GC offers insights into the shape of the Galactic potential (Kenyon et al. 2008).

HVSs within the Galactic bulge provide additional constraints on the characteristics of the stellar population in that particular region of the Galaxy. This includes the joint constraints of the stellar initial mass function of the GC and the HVS ejection rate, and the binary separation distribution (see Rossi, Kobayashi & Sari 2014; Rossi et al. 2017; Marchetti, Evans & Rossi 2022; Evans, Marchetti & Rossi 2022b).

Luna, Minniti & Alonso-García (2019) presented the first HVS candidates in the Galactic bulge using the Vista Variables in the Via Lactea (VVV) near-infrared (NIR) survey data. This work is extended here using the Gaia DR3 data in combination with a larger set of improved VVV-based NIR proper motions: VVV Infrared Astrometric Catalogue 2 (VIRAC2). We search for HVS candidates that can be further studied spectroscopically to confirm their nature and constrain their origin.

The paper is organized as follows: in Section 2, we describe the data sets. In Section 3, we present the selection of red clump stars and the computation of distances and tangential velocities. The resulting sample of high-velocity stars and HVSs candidates from the VIRAC2-Gaia DR3 crossmatch is presented in Section 4 and in Section 5 the equivalent sample of those stars having only VIRAC2 data. In Section 6, we summarize our conclusions.

2 THE Gaia DR3 AND VIRAC2 CATALOGUES

2.1 Gaia DR3

The Gaia DR3 provides, amongst other information, optical photometry (G, G_BP, and G_RP) and astrometry (parallaxes and proper motions) for ∼1.5 billion sources in the Galaxy, with 34 months of observations (Gaia Collaboration 2022). Gaia DR3 doubled the precision in proper motions and improved by 30 per cent the precision in parallaxes with respect to Gaia DR2. In the same vein, the errors improved by a factor of ∼2.5 in proper motions and up to 40 per cent in parallaxes (Lindegren et al. 2016; van Leeuwen et al. 2017; Gaia Collaboration 2021; Fabricius et al. 2021).

The survey has completeness below 60 per cent for sources at G ∼ 19 or fainter and in stellar densities of about |$5\times 10^5\, \mathrm{stars}\, \deg ^{-2}$|⁠, which is characteristic of crowded fields such as globular clusters. Its completeness further drops below 20 per cent even for brighter sources in fields with stellar densities of about |$10^6-10^7\, \mathrm{stars}\, \deg ^{-2}$| that are found in the Galactic bulge (Fabricius et al. 2021; Cantat-Gaudin et al. 2023). In the present study, we select Gaia DR3 sources that are matched with the VIRAC2 data (see next section) in coordinates.

2.2 VIRAC2

The VVV survey (Minniti et al. 2010) and its extension, the VVV eXtended survey (VVVX), acquired multi-epoch observations between 2010 and 2023 (hereafter we always refer to the survey as VVV, but it includes also the VVVX data). VVV produced a map of the Galactic bulge and southern part of the disc (⁠|$-130\deg \lt \, l\lt\, 20\deg$| and |$-15\deg \lt \, b\lt \, 10\deg$|⁠) covering |$\sim 1700\deg ^2$| in three NIR passbands: J(1.25 μm), H(1.64 μm), and K_s(2.14 μm) [the VVV original footprint is also covered in Z(0.87 μm) and Y(1.02 μm) bands]. For further information about the VVV survey and its data quality, we refer to publications by Minniti et al. (2010), Saito et al. (2012), Alonso-García et al. (2018), and Surot et al. (2019b). VIRAC2 (Smith et al., in preparation) is the second data release of the VIRAC (Smith et al. 2018). VIRAC2 is constructed from the PSF photometry catalogue, while VIRAC data are derived from aperture photometry. VIRAC2 is 90 per cent complete up to K_s ∼ 16 across the VVV bulge area (Sanders et al. 2022). Its proper motions are anchored to Gaia DR3 absolute reference frame. Luna et al. (2023) tested the reliability of VIRAC2 proper motions and their errors.

In our study, we use preliminary VIRAC2 data, and adopted the following quality selection criteria: (i) sources with a complete (5-parameter) astrometric solution, (ii) non-duplicates, (iii) sources that are detected in at least 20 per cent of the observations,² and (iv) sources with unit weight error (uwe) uwe < 1.2. The latter is used as a threshold to select single sources with good astrometric measurements.

3 SELECTION OF RED CLUMP SOURCES

Red clump (RC) stars are in the helium core-burning phase on their second ascent towards the giant branch (see Girardi 2016, for a review). They are easily identified in nearby stellar systems because they occupy a well-defined region in colour–magnitude diagrams (CMDs). They have a characteristic luminosity and exhibit minimal variability, which makes them a useful standard candle for distance determination (Alves et al. 2002; Ruiz-Dern et al. 2018).

Using the VIRAC2 data set, we select a region in the inner part of the Galactic bulge, in a box within |$-5\deg \lt \, l\lt \, 5\deg$| and |$-3\deg \lt \, b\lt \, 3\deg$|⁠. The total number of sources from VIRAC2 within the 60 deg² region is |$112\, 599\, 161$|⁠. In the same area, there are |$36\, 101\, 569$| sources in Gaia.

In Fig. 1, we show the J − K_s versus K_s CMD for all the stars in the 60 deg² bulge area centred on the GC. The blue sequence at J − K_s ≤ 0.6 mag, populated by the main-sequence disc, and the redder red giant branch (RGB) bulge sequence merge at K_s > 15 mag. The RGB is the nearly vertical sequence ranging in colour around J − K_s ∼ 0.8 − 1.1. Due to very high reddening in the centre-most regions of the bulge the RGB fans along the reddening vector covering the entire J − K_s range of the CMD. There are three overdensities along the RGB in the J − K_s CMD (Fig. 1). While the most obvious one, which we highlight inside an orange box in Fig. 1, corresponds to the RC stars from the bulge, the nature of the other two overdensities is more uncertain (Nataf et al. 2011, 2013; Alonso-García et al. 2018; Gonzalez et al. 2018).

In the CMD, we select the overdensity traced by the RC stars inside the orange box marked with a solid line in Fig. 1. These stars are selected within ±0.5 mag along the reddening vector

We have used the K_s total-to-selective ratio from (Alonso-García et al. 2017)

where the infrared colour excess of the RC stars is E(J − K_s) = (J − K_s) − (J − K_s)₀ and each star has an individual value taken from the reddening map of Sanders et al. (2022).

The orange box is centred on the RC. The intrinsic colour of the Galactic bulge RC is (J − K_s)₀ = 0.68 mag (Alves et al. 2002; Gonzalez et al. 2012; Ruiz-Dern et al. 2018), but the bulk of the RC in the bulge is at J − K_s ∼ 1.0 mag. The brightness of the RC can be determined from the K_s luminosity function (LF) shown in Fig. 2. The LF around the bulge RC is fitted with a Gaussian plus a polynomial (Gonzalez et al. 2013)

The fitted parameters are |$a=-4.33\times 10^5,~b=1.82\times 10^4,~c=2.6\times 10^6,~\sigma _{\mathrm{ RC}}=0.3\, \mathrm{ mag},~K_s^{\mathrm{ RC}}=13.27$| mag, where σ_RC = 0.3 mag results from the convolution of the intrinsic width of RC due to the distribution of stellar masses that populate the RC in the observed population (0.2 mag; Alves 2000), the observed width of the RC (0.22 mag; Surot et al. 2019a), and the photometric error from VIRAC2 (∼0.12 mag around the RC). In Fig. 2, the solid line is fit to the K_s LF following equation (3), and the dashed lines mark the limit of the RC selection. The LF for stars around (J − K_s) = 1 peaks at K_s = 13.3 mag. To exclude foreground main-sequence disc stars, we limit the RC selection polygon at the blue end to (J − K_s) ≥ 0.8 mag.

We then refine the selection of the RC candidates by a linear fit to the H − K_s versus J − K_s relation of the stars within the search box (⁠|$-5\deg \lt \, l\lt \, 5\deg$| and |$-3\deg \lt \, b\lt \, 3\deg$|⁠), and selecting ±0.25 mag to the sides of the fit. Fig. 3 shows the H − K_s versus J − K_s diagram, the polygon marks the boundary of the RC selection and the dashed line is the linear fit

After selecting the RC region in the CMD and in the colour–colour diagram from VIRAC2, we crossmatch the sources with Gaia DR3 with a 0|${_{.}^{\prime\prime}}$|5 tolerance in position, followed by removing spurious matches, which are identified as outliers in the distributions of (G − K_s) colours.

The mean separation of the crossmatched sources is 0.03 arcsec.

The Gaia DR3 – VIRAC2 matched sample consists of |$2\, 594\, 052$| RC candidate stars within the 60 deg² region in the inner bulge.

3.1 Foreground decontamination

The principal contaminant of the sample can be M-dwarf stars and reddened main-sequence disc stars and red giant stars from the bulge that populate the same area as the RC stars in the CMD and colour–colour diagrams (Alonso-García et al. 2018; Mejías et al. 2022). Those giant stars, at the faint end of the RC LF, belong to the RGB bump (Nataf et al. 2011), which overlaps with the more distant RC giants tracing the Milky Way structure behind the galactic bar (Gonzalez et al. 2018). Lim et al. (2021) studied the RC population in the southern bulge (⁠|$-6\,$| deg|$\lt \, b \, \lt -9\,$| deg) and found around 25 per cent of RGB bump interlopers in the magnitude range |$13\,$| mag|$\lt \, K_{s_0}\lt \, 14\,$| mag. A spectroscopy follow-up can discern the two populations (e.g. He, Luo & Chen 2022).

Bright foreground stars have well-characterized, reliable parallax in Gaia DR3 and/or in VIRAC2. Hence, we implement a cut in parallax removing those stars that lie within 5 kpc (ϖ > 0.2 mas) and are bright enough to have precise parallax (⁠|$\frac{\sigma _{\varpi }}{\varpi }\lt 0.5$|⁠). However, our sample might still contain faint foreground stars with poorly determined or no parallax measurements in the catalogues. We discuss this further in Section 4.1.

After the parallax cut our sample was reduced from |$2\, 594\, 052$| to |$2\, 077\, 026$| sources.

3.2 Population correction of RC absolute magnitude and distance determination

The absolute magnitude of RC depends on metallicity and age. The brightness of the RC in the Solar neighbourhood is |$M_{K_s}=-1.61$| (Alves et al. 2002; Ruiz-Dern et al. 2018; Sanders et al. 2022). To account for the metallicity and age dependence, a population correction is needed.

The metallicity distribution function (MDF) of RC stars in the bulge is bimodal (Zoccali et al. 2017; Queiroz et al. 2020, 2021; Rojas-Arriagada et al. 2020 ). The peaks are at around [Fe/H] = −0.4 and [Fe/H] = +0.3, with the contributions assumed to be 50 per cent each.³ We take from table 1 of Salaris & Girardi (2002) the stellar evolutionary model theoretical RC magnitudes for metallicity Z = 0.008 and Z = 0.03, corresponding to the two peaks, and average the contributions of these two bulge components to get the mean absolute K_s magnitude ‹|$M_{K_s}$|›.

We also consider three different ages: 8, 10, and 12 Gyr, which are typical of bulge stars (e.g. Surot et al. 2019a; Hasselquist et al. 2020; Sit & Ness 2020). For each age, we take the average ‹|$M_{K_{s_{Z1,Z2}}}$|› given the two metallicity components.

The difference between the local theoretical RC absolute magnitude – i.e. the absolute magnitude of the RC in the stellar evolutionary models of Salaris & Girardi (2002) – and the mean absolute RC magnitude for a given age and metallicity mix is the so-called population correction. It is derived from the Salaris & Girardi (2002) models in the following way:

where |$M_{K_s, \mathrm{ theo}}^{\mathrm{ RC}}=-1.54$|⁠.

The population correction (⁠|$\Delta M_{K_s}^{\mathrm{ RC}}$|⁠) variation between −0.011 and −0.099 mag, depending on age of the bulge RC stars, (see Table 1) is taken as the systematic error in the distance derivation.

The distance modulus μ₀ is then

where |$M_{K_s, \mathrm{ obs}.}^{\mathrm{ RC}}$| is the observationally determined absolute K_s magnitude for the local (Solar vicinity) RC population.

To compute the distance, we assume |$M_{K_s, \mathrm{ obs}.}^{\mathrm{ RC}}=-1.61\pm 0.07$| from Sanders et al. (2022), which is in agreement with Ruiz-Dern et al. (2018) |$M_{K_s}=(-1.606\pm 0.009)$|⁠. Then, applying the population correction described above, the distance in pc is expressed as

The heliocentric distances for our selection of RC candidate stars in the bulge CMD range between 7 and 11.5 kpc. For a typical RC star in our sample (K_s = 13.2 mag and extinction of |$A_{K_s}$| = 0.26 mag), the population correction implies a difference in distance of 130 pc compared to the distance derived only using the observed local absolute magnitude (⁠|$M_{K_s}=-1.54$| mag instead of |$M_{K_s}-1.61$| mag).

Adding to that, by assuming different ages, i.e. different corrections, the resulting distances vary by 250 pc, which at the distance of the bulge would be a relative error of 3 per cent. Such value is added to the error budget of the distance distribution.

3.3 Tangential velocity derivation

The high and variable dust extinction towards the Galactic bulge, added to the crowded stellar environment, might cause Gaia DR3 proper motions not be well-characterized (Lindegren et al. 2018; Riello et al. 2021; Fabricius et al. 2021; Battaglia et al. 2022; Luna et al. 2023). For those reasons, we make a separate selection based only on VIRAC2 proper motions, derived from NIR observations, for the following analysis.

After removing the bright foreground objects and having the distance estimation for the RC sources, the VIRAC2 proper motions (⁠|$\mu _{\alpha ^{*}}=\mu _{\alpha } \cos (\delta)$| and μ_δ) are corrected for the reflex motion of the Sun [(U_⊙, V_⊙, W_⊙) = (12.9, 245.6, 7.78) km s⁻¹] using the gala package (Price-Whelan 2017) from astropy (Astropy Collaboration 2013). For the correction computation we adopt the default Solar motion relative to the GC as a combination for the peculiar velocity (Schönrich, Binney & Dehnen 2010) and for the circular velocity at the Solar radius (Bovy 2015). The tangential velocity of each star in the sample is

where μ is the VIRAC2 proper motion and d is the heliocentric distance.

3.3.1 Monte Carlo sampling

We use a Monte Carlo sampling approach to derive a tangential velocity distribution for each star in our final sample. The VIRAC2 catalogue provides the correlation coefficient between the RA and Dec. components of the proper motions. Together with the proper motion uncertainties, one can use the covariance matrix of each source to generate a multivariate Gaussian distribution of the proper motion errors.

With that, we draw 100 000 random samples of proper motions for each star, which follow a Gaussian distribution. Each of those proper motions is corrected for the reflex motion of the Sun as described in the previous section.

For the distance, we assume a Gaussian distribution with a width given by the following parameters added in quadrature: the photometric error from VIRAC2 (∼0.12 mag around the RC), the absolute magnitude calibration error (0.07 mag; Sanders et al. 2022) and the population correction (0.088 mag; Salaris & Girardi 2002), resulting in a typical distance error of |$\sim 450\,$|pc.

From such distribution, we draw 100 000 random samples of distances for each source. Finally, for each source, we obtain 100 000 samples of tangential velocity, with the median of the distribution being their characteristic tangential velocity and the uncertainties are computed from the 16th and 84th percentile.

4 RESULTS AND DISCUSSION

The Gaia DR3-VIRAC2 crossmatch and foreground decontamination of the sample resulted in 2 077 026 sources located in the RC region of the NIR-CMD. We further restrict the selection of RC candidate stars to a narrower location in the CMD, within 0.3 < (J − K_s)₀ < 1 and −0.1 < (H − K_s)₀ < 0.4. From those, 49 sources have a lower limit in tangential velocity (16th percentile) above the local escape velocity, with a relative error in their tangential velocities smaller than 30 per cent (⁠|$\frac{\sigma _{V_{\mathrm{ tan}}}}{V_{\mathrm{ tan}}}\, \lt \, 0.3$|⁠) and |$\sigma _{V_{\mathrm{ tan}}}$| the difference between the 16th and 84th percentiles. The local escape velocity depends on the assumed potential, and for the Galactic bulge, it ranges between ∼600 and ∼800 |$\mathrm{km\, s^{-1}}$|⁠, reaching the latter within |$1\,$|kpc of the GC (Rossi et al. 2017; Deason et al. 2019).

To be more inclusive, we adopt the lower limit in velocity as the escape velocity given by the MWPotential2014 Galactic potential from the galpy package (Bovy 2015), with the addition of the influence of the bar – using the DehnenBarPotential from galpy – and a Kepler potential for Sgr A*. The potential is in the low-mass end of the Galaxy mass estimates (Callingham et al. 2019), with a total mass within 60 kpc of |$4\pm 0.7\times 10^{11}\, {\rm M}_{\odot }$| (Xue et al. 2008). In this potential, the dark matter halo is described by a NFW Potential, the bulge is modelled as a power-law density profile with an exponent of −1.8 and a cut-off radius of 1.9 kpc and the disc is modelled with a Miyamoto Nagai potential. The selection of the potential implies that some of the selected candidates move fast, but might still be bound, depending on their exact locations.

Their tangential velocities range between 731 and 1938 |$\mathrm{km\, s^{-1}}$|⁠, and are located between 0.3 and 3.5 kpc from the GC (considered to be at |$R_0=8.122\,$|kpc; GRAVITY Collaboration 2018). The range of tangential velocities is in agreement with the predictions from Generozov (2020) and data from Brown (2015), Boubert et al. (2018), and Koposov et al. (2020).

4.1 Possible sample contaminants and selection of interesting objects

The current data do not have sufficient information to assign a spectral type to each star. A precise spectroscopic determination of the stellar parameters (log g and T_eff) and the metallicity will enable us to confirm the stars as RC giants, hence confirming their location in the Galactic bulge.

To refine the sample of RC stars, we identify sources that might be M-dwarfs given their colours. For this, we follow the colour cuts by Mejías et al. (2022), based on a spectroscopically confirmed sample of M-dwarfs by West et al. (2011)

Three sources fall into the selection cuts, placing them as possible M-dwarfs. Their properties are listed in Appendix  A. Removing the likely M-dwarfs, the selection results in a sample of 46 stars. The excess of stars in the faint-end of the LF (solid histogram in Fig. 4) is expected. It is a consequence of a larger volume for the more distant, thus fainter, sources in the sample.

Fig. 5 shows the location of the 46 high-velocity stars in the extinction-corrected NIR CMD.

Given the heliocentric distance d derived from the RC calibration, we derive the distance with respect to the GC (r_GC) using astropy and the following parameters: the International Celestial Reference System (ICRS) coordinates of the GC (RA, Dec.) = (266.4051, −28.936175) deg (Reid & Brunthaler 2004); the velocity of the Sun in Galactocentric Cartesian coordinates (U_⊙, V_⊙, W_⊙) = (12.9, 245.6, 7.78) |$\mathrm{km\, s^{-1}}$| (Reid & Brunthaler 2004; Drimmel & Poggio 2018; GRAVITY Collaboration 2018); the distance from the Sun to the Galactic mid-plane z_⊙ = 20.8 pc (Bennett & Bovy 2019), and the distance from the Sun to the GC d_⊙ = 8.122 kpc (GRAVITY Collaboration 2018).

Fig. 6 shows the tangential velocity of the stars as a function of the galactocentric distance (r_GC). The solid line indicates the escape velocity at a given galactocentric distance, assuming the MWPotential2014 potential from the galpy package (Bovy 2015). The adopted potential is only a reference that helped us to select the threshold in velocity. Different assumptions would change the escape velocity by tens of |$\mathrm{km\, s^{-1}}$|⁠. As a reference for a higher escape velocity profile, the plot also shows the escape velocity profile given by the Irrgang13III potential from the galpy package. This potential is the Model III of Irrgang et al. (2013), with a total Galaxy mass within 50 kpc of |$8.1\times 10^{11}\, {\rm M}_{\odot }$|⁠.

The position in the sky of the final sample is shown in Fig. 7 in galactic coordinates. The background is the density plot of the VIRAC2 sources in the search box, where a whiter colour indicates a higher density, and is equivalent to a reddening map. The arrows represent the VIRAC2 proper motion vectors along the galactic latitude (b) and longitude (l) components, colour-coded by their tangential velocity. Black lines indicate the Gaia DR3 proper motion vector transformed to galactic coordinates and corrected for the reflex motion of the Sun. We note that the Gaia DR3 proper motion of some sources matches with that of VIRAC2, particularly for the sources with slower velocities.

In Fig. 7, we note asymmetries in the spatial distribution and proper motions of the sample. In the initial sample of Gaia DR3 – VIRAC2 crossmatch, the stars are distributed evenly in galactic latitude and longitude. However, in the sample of 46 high-velocity stars, 31 (67 per cent) stars are located at positive galactic latitudes. Furthermore, 39 (85 per cent) stars have VIRAC2 proper motion vectors pointing towards negative galactic longitude.

Based on the Gaia DR3 mock catalogue (Rybizki et al. 2020), a synthetic Milky Way catalogue that simulates the Gaia EDR3⁴ content, within the magnitude and colour range of our sample (17.8 < G < 19.7 mag; 2.6 < BP − RP < 3.9 mag), it is expected that ∼80 per cent of stars lie at a distance between 6.5 and 9.5 kpc. However, some of the stars in our final sample might still be foreground contaminants. Fig. 8 shows the tangential velocity of the stars as a function of their VIRAC2 proper motions. The dashed lines indicate the tangential velocity as a function of the proper motion for different fixed heliocentric distances, while the solid line is at the distance of the GC (8.122 kpc). Even if some of the stars were located at smaller heliocentric distances than our estimates, more than half of them are interesting objects travelling faster than the average disc rotation (246 |$\mathrm{km\, s^{-1}}$|⁠).

4.2 HVSs candidates and flight time

To test whether the HVS candidates radiate from the GC, we project the proper motion vectors, appropriately account for the Solar motion, and search for those that intersect a region within |$0.5\deg$| (∼75 pc at the distance of the GC) centred on the GC. This criterion is based on Speedystar⁵ (Evans, Marchetti & Rossi 2022b), where we simulate ejections of stars from the GC via the Hills mechanism and analyse the position angle of their proper motions corrected by the reflex motion of the Sun. The stars that are within 4 kpc from the GC, such as the stars in our sample, have a 2D velocity vector (projected proper motion) that approaches the GC in less than 0.45 deg. The criterion also accounts for the uncertainty in the GC distance (31 pc ∼ 0.2 deg), the GC position (0.66 mas; GRAVITY Collaboration 2018, 2021), and encloses the nuclear star cluster (0.03 deg; Schödel et al. 2014, 2023).

There are four sources that fulfil this condition. Their projected 2D velocity vector approaches the GC in a range between |$0.004$| and |$0.47\, \deg$|⁠; since these measurements lack radial velocity, the approach is not a physical quantity. The left panel of Fig. 9 shows their location in the CMD colour-coded by their tangential velocities. The top-right panel shows their tangential velocity as a function of their galactocentric distance colour-coded by their VIRAC2 proper motions, while the bottom-right panel shows their location in galactic coordinates, with the VIRAC2 vectors represented as arrows colour-coded by the tangential velocity.

To compute how long ago the HVS were ejected from the GC – the flight time, we project the VIRAC2 proper motion vectors backwards, in the direction towards the GC. We assume that the proper motion direction is not strongly affected by the projection in the sky

where t_f is the flight time in years, l and b are in degrees and |$\mu =\sqrt{\mu _{l^{*}}^2+\mu _b^2}$| in mas yr⁻¹.

The flight time of these four stars ranges from 1.1 × 10⁶ to 3.2 × 10⁶ yr. We can roughly estimate the ejection rate of RC stars by dividing the number of stars over the range of flight times, this is N_HVS/(max(t_F) − min(t_F)), where N_HVS is the observed number of RC HVS candidates, and the denominator is the difference between the maximum and minimum flight time of the candidates. This corresponds to an ejection rate of 1.9 × 10⁻⁶ yr⁻¹. However, we note that this is not a complete sample, and it is limited to the RC mass range, typically between 0.8 and 2 M_⊙ (Girardi 2016). Our sample is also affected by contamination from non-RC stars that can only be disentangled with a spectroscopic follow-up. Therefore, the ejection rates are lower than the total integrated rate, which current studies estimate at |$10^{-4}\, \mathrm{ yr}^{-1}$|⁠. This ejection rate estimate depends on the adopted IMF; if the IMF is top-light, the ejection rate can be |$10^{-2}\, \mathrm{ yr}^{-1}$|⁠, but if it is top-heavy, it can be |$10^{-4.5}\, \mathrm{ yr}^{-1}$| (e.g. Rossi et al. 2017; Marchetti, Evans & Rossi 2022; Evans, Marchetti & Rossi 2022a, b).

Table 2 shows for the four stars, the location in equatorial and galactic coordinates, and photometric information (K_s and J − K_s), Table 3 shows their VIRAC2 proper motions, and the derived parameters: tangential velocity (V_tan), galactocentric distance (r_GC), flight time (t_f), closest approach of the 2D velocity vector to the GC, and the probability (P) of exceeding the escape velocity, given by |$V_{\mathrm{ tan}, 16{\rm th}}(r_{\mathrm{ GC}})/V_{\mathrm{ esc}}(r_{\mathrm{ GC}})$|⁠, with |$V_{\mathrm{ tan}, 16{\rm th}}$| the lower limit in tangential velocity (16th percentile), and V_esc(r_GC) given by the galpyIrrgang13III potential; the stars that exceed such escape velocity have a probability of 1.

5 THE RED SAMPLE

In this section, we explore those sources that appear in VIRAC2, but not in Gaia DR3, because they are too faint or highly reddened.

We follow the procedure described in Sections 3 and 4, to select the RC stars, and clean the sample of foreground stars including M-dwarfs using VIRAC2 parallax and colour cuts.

The initial sample consists of 8 596 826 sources that are in the VIRAC2 catalogue, but do not have a crossmatch in Gaia DR3. From those, we remove stars closer than 5 kpc to us (ϖ > 0.2) and precise VIRAC2 parallax (⁠|$\frac{\sigma _{\varpi }}{\varpi }\lt 0.5$|⁠). The remaining sample has 8 214 069 sources. We select sources with lower limits in tangential velocity larger than the escape velocity at a given galactocentric distance. The resulting sample has 504 sources from which 481 have a relative error in their tangential velocities smaller than 30 per cent. From those sources, nine are possible M-dwarfs, identified using the colour cuts of equation (9). Thus, the final sample has 472 sources with tangential velocities up to V_tan = 2468 |$\mathrm{km\, s^{-1}}$|⁠. These stars are shown in Fig. 10. The left panel shows the location of the sources in the CMD corrected for extinction. The top-right panel shows the tangential velocity as a function of the galactocentric distance (r_GC). The bottom-right panel shows the location of the sources in galactic coordinates, colour-coded by their tangential velocity, and with the proper motion vector shown as arrows.

5.1 HVSs candidates and flight time

Fig. 11 is equivalent to Fig. 10, but shows the stars whose VIRAC2 proper motions point away from the GC. We note that in contrast to Fig. 7, the density of sources is higher towards lower galactic latitudes, where Gaia DR3 has fewer sources. Similar to Fig. 7, 72 per cent of the sources have VIRAC2 proper motions pointing towards negative galactic longitude, but they are distributed evenly in latitude with respect to the Galactic plane.

There are 65 sources that appear to be ejected from a region within |$0.5\deg$| from the GC, with proper motion vectors pointing out of it and with tangential velocities up to V_tan = 1875 |$\mathrm{km\, s^{-1}}$|⁠. Fig. 11 shows the 65 stars.

Their projected 2D velocity vector approaches the GC in a range between |$0.0004$| and |$0.49\, \deg$|⁠; as mentioned in Section 4.2, these values are not a physical quantity.

Following equation (10), the flight time of the 65 stars appearing to be ejected from the centre ranges between |$2.4\times 10^4$| and |$4.8\times 10^6\, \mathrm{ yr}$|⁠. This provides a rough estimate of an ejection rate of |$1.4\times 10^{-5}\, \mathrm{ yr}^{-1}$|⁠, an order of magnitude higher than for the sample that includes also Gaia DR3 measurements. As explained in Section 4.2, this estimate is lower than the integrated ejection rate. Fig. 12 shows the distribution of flight times, with the major ejection episode occurring before 2 Myr ago.

Table 4 is an extract of Table B1 and shows for the 65 stars with VIRAC2 and no Gaia DR3 counterpart, the location in equatorial and galactic coordinates, and photometric information (K_s and J − K_s). Table 5 is an extract of Table B2 and shows VIRAC2 proper motions, and the derived parameters: tangential velocity (V_tan), galactocentric distance (r_GC), flight time (t_f), the closest approach of the 2D velocity vector to the GC, and the probability (P) of exceeding the escape velocity.

HVSs can also originate from a globular cluster hosting an intermediate-mass black hole (IMBH). Stellar interactions with the IMBH might result in the ejection of HVSs in a similar way as the Hills mechanism (e.g. Fragione & Gualandris 2019). In that case, the star would point back to its parent cluster.

The spatial and velocity distributions can be used to discern those stars from (hyper)runaway stars originating from the Galactic disc, where the runaway stars have lower velocities and are closer to the Galactic plane (e.g. Fragione & Gualandris 2019; Generozov & Perets 2022). In the same way, the velocity distribution of HVSs originating from the GC would have a larger contribution to high-velocity stars.

A bMBH, like a pair of SMBH–IMBH, would produce HVSs with potentially anisotropic velocity distribution, depending on the eccentricity and mass ratio of the binary stellar system that gets disrupted (e.g. Quinlan 1996; Dotti et al. 2006; Haardt, Sesana & Madau 2006; Sesana, Haardt & Madau 2006). The velocity and spatial distribution predicted by the Hills mechanism are isotropic, thus, the anisotropy in the velocity and spatial distribution of HVSs originating from a bMBH can potentially distinguish their ejection scenario from the SMBH interaction.

The spatial distribution predicted by the Hills mechanism is isotropic. In Fig. 13, we show the position angle (PA) of the sample of 65 stars that point away from the GC. There is a peak in the cos(PA) since most of the stars point towards negative longitudes, but without a trend with velocity (see Fig. 14). Hence, the direction is not isotropic, but the velocity is.

Fig. 15 shows the inverse cumulative distribution function of the high-velocity stars’ tangential velocity in both samples: black is the distribution of sources matched in VIRAC2 and Gaia DR3 and that for sources without Gaia DR3 counterpart is plotted with red line. The stars in the sample with only VIRAC2 data are closer to the GC. This sample also contains stars that reach higher velocities than the sample with Gaia DR3 counterpart. With the deceleration occurring in the first |$\sim 200\,$|pc (Kenyon et al. 2008), the stars with Gaia DR3 counterpart, that lie at larger distances from the GC, might have been decelerated.

6 CONCLUSIONS

We have identified 69 HVS candidates in the Milky Way bulge that appear to be ejected from the GC. From those, four have a Gaia DR3 counterpart, but with Gaia proper motions that are not reliable and 65 were detected only in the VIRAC2 sample. All these stars are selected among the likely RC candidates that place them within the bulge. Their flight time ranges between 2.4 × 10⁴ and 4.8 × 10⁶ yr, and implies an integrated ejection rate of the order of 1.4 × 10⁻⁵ yr⁻¹, provided that they all originated from the GC and are currently still within the bulge. Nevertheless, since the surveys are not complete, our sample is not complete. Adding to that, it is limited to RC stars, thus the derived ejection rates are lower than the current estimates (⁠|$10^{-4}\, \mathrm{ yr}^{-1}$|⁠).

To fully determine the HVS candidates’ spatial movement and reconstruct their orbits, we need to obtain the missing radial velocities. The spectroscopic follow-up will also provide the chemical composition and spectral type measurements necessary to confirm the candidates as RC stars. In the case a given star is not an RC star, the distance estimation would not be correct, nor its tangential velocity. However, such stars might be still interesting as they would be moving faster than the average disc rotation. Hence, they would offer the possibility of exploring the production scenarios of such stars through e.g. binary interactions or supernova kicks.

Having the first detection of HVSs in the Galactic bulge confirmed by future spectroscopy, we will be in a position to constrain the existing models, compare the observed with HVSs mock population (Marchetti et al. 2018) and compute their orbits (e.g. Fernández-Trincado et al. 2020; GravPot16⁶), and thus probe if they come from the GC or have another origin.

The Hills mechanism is parametrized by the shape of the stellar initial mass function in the GC, the distribution of binary orbital periods and mass fraction, and the ejection rate. The number of confirmed HVSs and their properties will constrain primarily the ejection rate and the shape of the IMF in the adopted model (e.g. Evans, Marchetti & Rossi 2022b, Speedystar). Interestingly, the discovery of the first HVSs close to the production centre will give a unique possibility to explore a recent stellar interaction with Sgr A*.

The velocity distribution is compatible with current ejection models, where the high-velocity stars should reside close to the GC, and those ejected with the Hills mechanism will also have an isotropic spatial distribution. A different spatial distribution, or proper motions that do not point back to the GC can be explained by other ejection scenarios, such as a disruption of globular clusters by the SMBH, inspiralling IMBH, dynamic interactions with binary massive BH, or supernova kicks. The high-velocity stars point back to the production origin, being these globular clusters, satellite galaxies such as the Magellanic Clouds, or the Galactic disc, where the runaway scenario is more probable.

The future Gaia data releases⁷ will improve the astrometry precision by a factor of ∼2.5, reaching 0.5 and 0.3 mas yr⁻¹ in DR4 and DR5, respectively, for the faintest stars (G = 20.7 mag). The accuracy will improve as well, as the baseline of observations will increase to 60 months in DR4 and 120 months in DR5. The radial velocities are limited to sources brighter than G ∼ 15.5 mag and the spectra at G ∼ 17.5 mag. Thus, to confirm the nature of the candidate HVSs in the Galactic bulge, a dedicated NIR spectroscopic follow-up is still necessary. GaiaNIR is a proposed space mission that will map the entire sky in the NIR, extending the capabilities and deepness of Gaia; if accepted, it will be launched in the 2040s, extending the baseline by 20 yr, which translates into a 20 × more accurate proper motions. The future ground and space facilities might change that scenario. The Extremely Large Telescope (ELT) will reach μas yr⁻¹ precision in proper motions (Trippe et al. 2010; Rodeghiero et al. 2021), and there are JWST observations focusing on the Galactic bulge. The Wide-Field Instrument onboard the Roman telescope will be ideal for NIR surveys with precise photometry (PSF full width at half-maximum <0.1 arcsec) and low resolution (R ∼ 600) multi-object spectroscopy, sufficient for radial velocity measurements. The Rubin observatory will map the southern sky and will reach proper motion precision of 0.2(1) mas yr⁻¹ for sources with magnitude r = 21(24) mag.⁸

ACKNOWLEDGEMENTS

AL acknowledges support from the Agencia Nacional de Investigación y Desarrollo (ANID) Doctorado Nacional 2021 scholarship 21211520, and the ESO studentship. NWCL gratefully acknowledges the generous support of a Fondo Nacional de Desarrollo Científico y Tecnológico (FONDECYT) regular grant number1230082, as well as support from Millenium Nucleus NCN19_058 (TITANs) and funding via the BASAL Centro de Excelencia en Astrofisica y Tecnologias Afines (CATA) grant PFB-06/2007. NWCL also thanks for support from ANID BASAL project number ACE210002 and ANID BASAL project numbers ACE210002 and FB210003. JA-G acknowledges support from FONDECYT regular grant number 1201490 and ANID – Millennium Science Initiative Program – ICN12_009 awarded to the Millennium Institute of Astrophysics MAS. AVN acknowledges support from the ANID scholarship programme Doctorado Nacional 2020-21201226, ANID, Millennium Science Initiative, ICN12_009 and ANID BASAL FB210003. DM gratefully acknowledges support by the ANID BASAL project numbers ACE210002 and FB210003, by FONDECYT project number 1220724, and by CNPq/Brazil through project number 350104/2022-0. We gratefully acknowledge the use of data from the ESO Public Survey programme IDs 179.B-2002 and 198.B-2004 taken with the VISTA telescope and data products from the Cambridge Astronomical Survey Unit. 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 2013, 2018). This research made use of numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), and scikit-learn (Pedregosa et al. 2011).

DATA AVAILABILITY

We use data from the ESO Public Survey programme IDs 179.B-2002 and 198.B-2004 taken with the VISTA telescope and data products from the Cambridge Astronomical Survey Unit. The VVV reduced images are available in the ESO Science Archive. 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). For the HVS candidates, we provide tables in the paper that will be available in CDS.

Footnotes

References

APPENDIX A: POSSIBLE M-DWARFS IN THE Gaia DR3 – VIRAC2 CROSSMATCH

In this section, we present the list of the possible M-dwarfs in our Gaia DR3 – VIRAC2 sample (Table A1) and in the VIRAC2 sample without Gaia DR3 counterparts (Table A2). These stars are selected with the colour cuts of equation (9), as described in Section 4.1.

APPENDIX B: STARS IN VIRAC2 POINTING AWAY FROM THE GC

This section presents the photometric properties and derived parameters of the 65 stars that have VIRAC2 data, but no Gaia DR3 counterpart, and whose VIRAC2 proper motions point away from |$0.5\deg$| from the GC.