New candidate hypervelocity red clump stars in the inner Galactic bulge

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 NIR 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 4 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.


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 (Barbuy et al. 2018;Springel & Hernquist 2005;Steinmetz & Muller 1995).The environment in such region is extreme, with high stellar crowding and in its centre resides a supermassive black hole (SMBH), Sgr A*, with a mass of 4.3 × 10 6 M ⊙ (GRAVITY Collaboration et al. 2022;Genzel et al. 2010;Ghez et al. 2008;Rodriguez 1978).The high stellar density favours interactions amongst stars and between stars and SgrA*.
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 ★ E-mail: alonso.luna@eso.orgdynamics surrounding the interaction itself.For example, Blaauw (1961) showed that runaway stars with peculiar velocities larger than 30 km s −1 can form after a binary stellar system disruption in the supernova explosion of one of its components.Poveda et al. (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 km s −1 and in some cases reaching nearly 200 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 SgrA*; one star is ejected as an HVS, while the other remains attached to the SMBH as an S-Star 1 (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 a SMBH (Brown 2015;Capuzzo-Dolcetta & Fragione 2015;Fragione et al. 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:    ∝ () −1/2 (/ ⊙ ) 1/6 .Given the mass of Sgr A*, 4 × 10 6 M ⊙ , HVSs can be ejected at up to ∼ 4 000 km s −1 (Hills 1988;Rossi et al. 2017Rossi et al. , 2021;;Sari et al. 2010), which is larger than the escape velocity.The local escape velocity is ∼ 830 km s −1 in the GC, ∼320 km s −1 in the halo (50 kpc) and ∼530 km s −1 at the Sun location (Deason et al. 2019;Rossi et al. 2017).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* (Generozov 2021;Hills 1988).
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 km s −1 , exceeding the local escape speed (e.g. Brown et al. 2018;Kenyon et al. 2014;Kollmeier et al. 2009;Palladino et al. 2014).Almost all of them are B-type stars with masses between 2.5 and 4 ⊙ , except for the fastest star amongst them: HVS-S5 (Koposov et al. 2020), an A-type star with a velocity of 1700 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.Boubert et al. 2020;Evans et al. 2021;Lennon et al. 2017;Przybilla et al. 2008), 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 dwarfs, travelling at space velocities larger than 1 300 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 Evans et al. 2022b;Marchetti et al. 2022;Rossi et al. 2014Rossi et al. , 2017)).Luna et al. (2019) presented the first HVS candidates in the Galactic bulge using the VVV near IR 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: VIRAC2.We search for HVS candidates that can be further studied spectroscopically to confirm their nature and constrain their origin.
The paper is organised 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 summarise our conclusions.

Gaia DR3
The Gaia third data release (DR3) provides, amongst other information, optical photometry (,   and   ) and astrometry (parallaxes and proper motions) for ∼1.5 billion sources in the Galaxy, with 34 months of observations (Gaia Collaboration et al. 2022).Gaia DR3 doubled the precision in proper motions and improved by 30% 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% in parallaxes (Fabricius et al. 2021a;Gaia Collaboration et al. 2021;Lindegren et al. 2016;van Leeuwen et al. 2017).
The survey has completeness below 60% for sources at  ∼ 19 or fainter and in stellar densities of about 5 × 10 5 stars deg −2 , which is characteristic of crowded fields such as globular clusters.Its completeness further drops below 20% even for brighter sources in fields with stellar densities of about 10 6 − 10 7 stars deg −2 that are found in the Galactic bulge (Cantat-Gaudin et al. 2023;Fabricius et al. 2021b).In the present study, we select Gaia DR3 sources that are matched with the VIRAC2 data (see next section) in coordinates.

VIRAC2
The VVV survey (Minniti et al. 2010) and its extension, the VVV eXtended survey (VVVX), acquired multi-epoch observations between 2010 to 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 <  < 20 deg and −15 deg <  < 10 deg) covering ∼ 1700 deg 2 in three near-infrared passbands:  (1.25 ),  (1.64 ) and   (2.14 ) (the VVV original footprint is also cover in  (0.87 ) and  (1.02 ) 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), andSurot et al. (2019b).VIRAC2 (Smith et al., in prep.) is the second data release of the VVV Infrared Astrometric Catalogue (VIRAC) (Smith et al. 2018).VIRAC2 is constructed from the PSF photometry catalogue, while VIRAC data are derived from aperture photometry.VIRAC2 is 90% complete up to   ∼ 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 (5parameter) astrometric solution, (ii) non-duplicates, (iii) sources that are detected in at least 20% of the observations, 2 , and (iv) sources with unit weight error (uwe)  < 1.2.The latter is used as a threshold to select single sources with good astrometric measurements.

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 <  < 5 deg and −3 deg <  < 3 deg.The total number of sources from VIRAC2 within the 60 sq.deg region is 112 599 161.In the same area, there are 36 101 569 sources in Gaia .
In Figure 1 we show the  −   vs.   CMD for all the stars in the 60 sq.deg bulge area centred on the GC.The blue sequence at  −   ≤ 0.6 mag, populated by the Main Sequence disc, and the redder Red Giant Branch (RGB) bulge sequence merge at   > 15 mag.The RGB is the nearly vertical sequence ranging in colour around  −   ∼ 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  −   range of the CMD.There are three overdensities along the Red Giant Branch (RGB) in the  −   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 (Alonso-García et al. 2018;Gonzalez et al. 2018;Nataf et al. 2011Nataf et al. , 2013)).
In the CMD, we select the overdensity traced by the RC stars inside the orange box marked with a solid line in Figure 1.These stars are selected within ±0.5mag along the reddening vector:   = 12.7 − 0.428( −   ) (1) We have used the   total-to-selective ratio from (Alonso-García et al. 2017): where the infrared colour excess of the RC stars is  ( −   ) = ( −   ) − ( −   ) 0 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 ( −   ) 0 = 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  −   ∼ 1.0 mag.The brightness of the RC can be determined from the   luminosity function (LF) shown in Figure 2. The LF around the bulge RC is fitted with a Gaussian plus a polynomial (Gonzalez et al. 2013): The fitted parameters are:  = −4.33 × 10 5 ,  = 1.82 × 10 4 ,  = 2.6 × 10 6 ,   = 0.3,    = 13.27mag.; where   = 0.3 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 Figure 2, the solid line is fit to the   LF following Equation 3, and the dashed lines mark the limit of the RC selection.The LF for stars around ( −   ) = 1 peaks at   = 13.3 mag.To exclude foreground Main Sequence disc stars, we limit the RC selection polygon at the blue end to ( −   ) ≥ 0.8 mag.
We then refine the selection of the RC candidates by a linear 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. ′′ 5 tolerance in position, followed by removing spurious matches, which are identified as outliers in the distributions of ( −   ) 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 sq.deg region in the inner bulge.

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 colourcolour 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<  < −9 deg) and found around 25% of RGB bump interlopers in the magnitude range 13 mag<   0 < 14 mag.A spectroscopy follow-up can discern the two populations (e.g., He et al. 2022).
Bright foreground stars have well-characterised, 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 (    < 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.

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    = −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 (Queiroz et al. 2020(Queiroz et al. , 2021;;Rojas-Arriagada et al. 2020;Zoccali et al. 2017).The peaks are at around [/] = −0.4 and [/] = +0.3, with the contributions assumed to be 50% each3 .We take from Table 1 of Salaris & Girardi (2002) the stellar evolutionary model theoretical RC magnitudes for metallicity  = 0.008 and  = 0.03, corresponding to the two peaks, and average the contributions of these two bulge components to get the mean absolute   magnitude <    >.
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     , ℎ = −1.54.

The population correction (Δ𝑀 𝑅𝐶 𝐾 𝑠
) 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  0 is then: where     ,. is the observationally determined absolute   magnitude for the local (Solar vicinity) RC population.
To compute the distance, we assume     ,.= −1.61± 0.07 from Sanders et al. (2022), which is in agreement with Ruiz-Dern et al. ( 2018)    = (−1.606±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 (  = 13.2 mag and extinction of    =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 (   = −1.54mag instead of    − 1.61 mag).
Adding to that, by assuming different ages, that is, different corrections, the resulting distances vary by 250 pc, which at the distance of the bulge would be a relative error of 3%.Such value is added to the error budget of the distance distribution.

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 characterised (Battaglia et al. 2022;Fabricius et al. 2021b;Lindegren et al. 2018;Luna et al. 2023;Riello et al. 2021).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 (  * =   cos() and   ) are corrected for the reflex motion of the Sun (( ⊙ ,  ⊙ ,  ⊙ ) = (12.9,245.6, 7.78) km s −1 ) using the gala package (Price-Whelan 2017) from Astropy (Astropy Collaboration et al. 2013a).For the correction computation we adopt the default Solar motion relative to the GC as a combination for the peculiar velocity (Schönrich et al. 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  is the heliocentric distance.

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 ∼ 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 16 th and 84 th percentile.

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 < ( −   ) 0 < 1 and −0.1 < ( −   ) 0 < 0.4.From those, 49 sources have a lower limit in tangential velocity (16 th percentile) above the local escape velocity, with a relative error in their tangential velocities smaller than 30% (        < 0.3) and     the difference between the 16 th and 84 th percentiles.The local escape velocity depends on the assumed potential, and for the Galactic bulge, it ranges between ∼ 600 and ∼ 800 km s −1 , reaching the latter within 1 kpc of the GC (Deason et al. 2019;Rossi et al. 2017).
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 SgrA*.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 ± 0.7 × 10 11  ⊙ (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 km s −1 , and are located between 0.3 and 3.5 kpc from the GC (considered to be at  0 = 8.122 kpc (GRAVITY Collaboration et al. 2018)).The range of tangential velocities is in agreement with the predictions from Generozov (2020) and data from Boubert et al. (2018); Brown (2015); Koposov et al. (2020).

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  and  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)   Figure 5 shows the location of the 46 high-velocity stars in the extinction-corrected NIR CMD.
Figure 6 shows the tangential velocity of the stars as a function of the galactocentric distance (  ).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 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 × 10 11  ⊙ .
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 () and longitude () 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%) stars are located at positive galactic latitudes.Furthermore, 39 (85%) 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 EDR34 content, within the magnitude and colour range of our sample (17.8 <  < 19.7 mag; 2.6 <  −  < 3.9 mag), it is expected that ∼80% 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.Figure 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 km s −1 ).

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  Speedystar5 (Evans et al. 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.2deg), the GC position (0.66 mas) (GRAVITY Collaboration et al. 2018Collaboration et al. , 2021)), and encloses the nuclear star cluster (0.03 deg) (Schödel et al. 2014(Schödel et al. , 2023)).
There are 4 sources that fulfil this condition.Their projected 2D velocity vector approaches the GC in a range between 0.004 deg 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 colourcoded 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   is the flight time in years,  and  are in degrees and in mas yr −1 .The flight time of these 4 stars ranges from 1.1 × 10 6 yr to 3.2 × 10 6 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    /((  ) − (  )), where    is the observed number of Red Clump 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 −6 yr −1 .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  ⊙ (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  −1 .This ejection rate estimate depends on the adopted IMF; if the IMF is top-light, the ejection rate can be 10 −2  −1 , but if it is top-heavy, it can be 10 −4.5  −1 (e.g., Evans et al. 2022a,b;Marchetti et al. 2022;Rossi et al. 2017).
Table 2 shows for the 4 stars, the location in equatorial and galactic coordinates, and photometric information (  and  −   ), Table 3 shows their VIRAC2 proper motions, and the derived parameters: tangential velocity (   ), galactocentric distance (  ), flight time (  ), closest approach of the 2D velocity vector to the GC, and the probability () of exceeding the escape velocity, given by   ,16 th (  )/  (  ), with   ,16 th the lower limit in tangential velocity (16 th percentile), and   (  ) given by the galpy Irrgang13III potential; the stars that exceed such escape velocity have a probability of 1.

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 Sec. 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 (    < 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%.From those sources, 9 are possible M-dwarfs, identified using the colour cuts of Eq.9.Thus, the final sample has 472 sources with tangential velocities up to    = 2468 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 (  ).The bottom-right panel shows the location of the sources in galactic coordinates, colourcoded by their tangential velocity, and with the proper motion vector shown as arrows.

HVSs candidates and flight time
Figure 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% 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    = 1875 km s −1 .Fig. 11 shows the 65 stars.
Their projected 2D velocity vector approaches the GC in a range between 0.0004 deg and 0.49 deg; as mentioned in Section 4.2, these values are not a physical quantity.
Following Eq.10, the flight time of the 65 stars appearing to be ejected from the centre ranges between 2.4 × 10 4  and 4.8 × 10 6 .This provides a rough estimate of an ejection rate of 1.4 × 10 −5  −1 , an order of magnitude higher than for the sample that includes also Gaia DR3 measurements.As explained in Sec.4.2, this estimate is lower than the integrated ejection rate.Figure 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 (  and  −   ).Table 5 is an extract of Table B2 and shows VIRAC2 proper motions, and the derived parameters: tangential velocity (   ), galactocentric distance (  ), flight time (  ), the closest approach of the 2D velocity vector to the GC, and the probability () 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 highvelocity stars.
A binary massive black hole (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., Dotti et al. 2006;Haardt et al. 2006;Quinlan 1996;Sesana et al. 2006).The velocity and spatial Table 2. Location and photometric properties of the stars with Gaia DR3 counterpart and VIRAC2 proper motions pointing away from 0.5 deg around the GC.The Galactic longitude and latitude have a precision of 0. ′′ 5. 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 ( ) 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.
Figure 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 ∼ 200 pc (Kenyon et al. 2008), the stars with Gaia DR3 counterpart, that lie at larger distances from the GC, might have been decelerated.

CONCLUSIONS
We have identified 69 HVS candidates in the MW bulge that appear to be ejected from the GC.From those, 4 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 4 and 4.8 × 10 6 yr, and implies an integrated ejection rate of the order of 1.4 × 10 −5 yr −1 , 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  −1 ).
To fully determine the HVS candidates' spatial movement and reconstruct their orbits, we need to obtain the missing radial velocities.  5. Derived parameters of the stars with no Gaia DR3 counterpart and VIRAC2 proper motions pointing out from 0.5 deg around the GC.The second to last column indicates how close the projected 2D velocity vector approaches the GC.The last column is the probability of exceeding the escape velocity.This is an extract, the complete table can be found in Appendix B. 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, for example, 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 6 ), and thus probe if they come from the GC or have another origin.
The Hills mechanism is parameterised by the shape of the stellar 6 https://gravpot.utinam.cnrs.frinitial 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 et al. 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 7 will improve the astrometry precision by a factor of ∼ 2.5, reaching 0.5 mas yr −1 and 0.3 mas yr −1 in DR4 and DR5, respectively, for the faintest stars ( = 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  ∼ 15.5 mag and the spectra at  ∼ 17.5 mag.Thus, to confirm the nature of the candidate HVSs in the Galactic bulge, a dedicated near-IR 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 years, which translates into a 20× more accurate proper motions.The future ground and space facilities might change that scenario.The ELT will reach as yr −1 precision in proper motions (Rodeghiero et al. 2021;Trippe et al. 2010), and there are JWST observations focusing on the Galactic bulge.The Wide-Field Instrument on board the Roman telescope will be ideal for NIR surveys with precise photometry (PSF FWHM < 0. ′′ 1 ) and low resolution ( ∼ 600) multi-object spectroscopy, 7 https://www.cosmos.esa.int/web/gaia/science-performancesufficient for radial velocity measurements.The Rubin observatory will map the southern sky and will reach proper motion precision of 0.2(1) mas yr −1 for sources with magnitude  = 21(24) mag8 .

Figure 1 .
Figure 1.VVV near-IR CMD of the 60 sq.deg studied region towards the Galactic bulge.The orange box encloses the area used for the RC sources selection.

Figure 2 .
Figure 2.   luminosity function around the less reddened RC location of the RGB selected around (  −   ) = 1.The solid line is fit to the distribution and follows Equation 3, and the dashed lines are the adopted limits of the RC selection.

Figure 3 .
Figure 3. Colour-colour plot covering the 60 sq.deg area search box (−5 deg <  < 5 deg and −3 deg <  < 3 deg).The dashed line is the linear fit to refine the RC selection.The orange box encloses further selection of RC stars.
, based on a spectroscopically confirmed sample of M dwarfs by West et al. (2011): 0.414 <  −  < 0.695 0.058 <  −   < 0.504 0.621 <  −   < 1.051 (9) 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

Figure 4 .
Figure 4.   0 luminosity function.The y-axis shows the normalised count of sources.The dashed histogram are the stars within the magnitude range 12.5 <   0 < 13.75.The solid histogram are the 46 stars that travel with a tangential velocity larger than the local escape velocity and that are not identified as M-dwarfs by their colour.

Figure 5 .
Figure 5. Deredenned CMD.The colored points are stars that travel with a tangential velocity larger than the local escape velocity.

Figure 6 .
Figure 6.Tangential velocity vs galactocentric distance (  ).The solid and dashed lines are the escape velocity at a given   assuming the galpy Galactic potential MWPotential2014, or Irrgang13III, respectively.. The coloured points are the 265 stars remaining after foreground objects' decontamination, and colour and quality cuts.The points are colour-coded by their VIRAC2 proper motions (PM).

Figure 7 .
Figure 7. Stars that travel with a tangential velocity larger than the local escape velocity.The arrows represent the VIRAC2 proper motion vectors along the galactic latitude () and longitude () components after correcting for the reflex motion of the Sun.The points and arrows are colour-coded by the tangential velocity of the source.Black lines indicate the Gaia DR3 proper motion vector transformed to galactic coordinates and corrected for the reflex motion of the Sun.The cross places the GC at (, ) = (0, 0) and the circle surrounding it represents a radius of 0.5 deg ≃ 75 pc.

Figure 8 .
Figure 8. Tangential velocity as a function of proper motion at different fixed heliocentric distances.The blue points are the 46 high-velocity stars in the sample described in Sect.4.1, their proper motions are from VIRAC2 and the tangential velocities are computed with Eq. 8.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.

Figure 9 .
Figure 9. Near-infrared photometry, position and kinematics of the candidate HVSs detected in both VIRAC2 and Gaia DR3.Stars that appear to be ejected from within 0.5 deg ≃ 75 pc of the GC.The left panel is the CMD, the same as in Fig. 5.The top-right panel is the tangential velocity as a function of galactocentric distance, same as Fig.6.The bottom-right panel is the location of the stars in galactic coordinates with the VIRAC2 proper motions shown as arrows colour-coded by their tangential velocities, same as Fig. 7.In the bottom-right panel, the cross marks the GC, and the surrounding circle has a radius of 0.5 deg ≃ 75 pc.

Figure 10 .
Figure 10.Near-infrared photometry, position and kinematics of the RC stars detected in VIRAC2, but not in Gaia DR3.The left panel is the CMD, same as in Fig. 5; the top-right panel is the tangential velocity as a function of galactocentric distance, same as Fig.6; and the bottom-right panel is the location of the stars in galactic coordinates with the VIRAC2 proper motions shown as arrows colour-coded by their tangential velocities, same as Fig. 7.

Figure 11 .
Figure 11.Near-infrared photometry, position and kinematics of the candidate HVSs detected in VIRAC2, but not in Gaia DR3.Stars that appear to be ejected from within 0.5  ≃ 75  of the GC.The left panel is the CMD, same as in Fig. 5; the top-right panel is the tangential velocity as a function of galactocentric distance, same as Fig.6; and the bottom-right panel is the location of the stars in galactic coordinates with the VIRAC2 proper motions shown as arrows colour-coded by their tangential velocities, same as Fig. 7.

Figure 12 .
Figure 12.Flight time in Myr for the candidate HVSs that appear in VIRAC2, but not in Gaia DR3 and have a proper motion vector pointing away from a location within 1 deg ≃ 75 pc of the GC.

Figure 13 .
Figure 13.Cosine of the position angle of the stars in VIRAC2 that point out from 0.5 deg around the GC.

Figure 14 .
Figure 14.Tangential velocity as a function of the Cosine of the position angle of the stars in VIRAC2 that point out from 0.5 deg around the GC.

Figure 15 .
Figure 15.Inverse cumulative distribution of the high-velocity stars.The black histogram is for the final sample of stars in the VIRAC2 and Gaia DR3 crossmatch (solid) and those that point out from 0.5 deg around the GC (dashed).The black histogram is for the final sample of stars in VIRAC2 (solid) and those that point out from 0.5 deg around the GC (dashed).

Table 1 .
Population correction (Δ   ) for the   absolute magnitude of RC stars.Δ   is computed with Table 1 of Salaris & Girardi (2002) assuming different ages and metallicities of Galactic bulge stars.

Table 3 .
(Rybizki et al. 2022f the stars with Gaia DR3 counterpart and VIRAC2 proper motions pointing out from 0.5 deg around the GC.The second to last column indicates how close the projected 2D velocity vector approaches the GC.The last column indicates the probability of exceeding the escape velocity.*Accordingto the classifier of spurious astrometric solution in Gaia DR3(Rybizki et al. 2022), these sources have a good astrometric solution.The VIRAC2 proper motions are consistent within 1.

Table 4 .
Location and photometric properties of the stars with no Gaia DR3 counterpart and VIRAC2 proper motions pointing away from 0.5 deg around the GC.The Galactic longitude and latitude have a precision of 0. ′′ 5.This is an extract, the complete table can be found in Appendix B.

Table A1 .
Properties of the stars, in the Gaia DR3 -VIRAC2 crossmatch, that could be M-dwarfs according to the colour cuts of Eq. 9.The proper motions are from VIRAC2, without the reflex motion correction.

Table A2 .
Properties of the stars, in VIRAC2 without Gaia DR3 counterpart, that could be M-dwarfs according to the colour cuts of Eq. 9.The proper motions are from VIRAC2, without the reflex motion correction.  [mag]  −   [mag] ]