Improved constraints on Galactic Centre ejection of hypervelocity stars based on novel search method

ABSTRACT

Hypervelocity stars (HVSs) are stars which have been ejected from the Galactic Centre (GC) at velocities of up to a few thousand |$\text{km}\, \text{s}^{-1}$|⁠. They are tracers of the Galactic potential and can be used to infer properties of the GC, such as the initial mass function and assembly history. HVSs are rare, however, with only about a dozen promising candidates discovered so far. In this work, we make use of a novel, highly efficient method to identify new HVS candidates in Gaia. This method uses the nearly radial trajectories of HVSs to infer their distances and velocities based on their position and Gaia proper motion alone. Through comparison of inferred distances with Gaia parallaxes and photometry, we identified 600 HVS candidates with G < 20 including the previously discovered S5-HVS1, out of which we obtained ground-based follow-up observations for 196 stars. As we found no new HVSs based on their radial velocity, we used detailed HVS ejection simulations to significantly improve previous HVS ejection rate constraints. In particular, the ejection rate of HVSs more massive than 1 |$\mathrm{M_\odot }$| cannot be higher than |$10^{-5}$| yr|$^{-1}$| at |$2\sigma$| significance. Additionally, we predict that there are 5–45 unbound HVSs in the complete Gaia catalogue (⁠|$1\sigma$| interval), most of which will be main-sequence stars of a few M|$_\odot$| at heliocentric distances of tens to hundreds of kpc. By comparing our results to literature HVS candidates, we find an indication of either a time-dependent ejection rate of HVSs or a non-GC origin of previously identified HVS candidates.

1 INTRODUCTION

The Galactic Centre (GC) is a highly complex environment, hosting both a nuclear star cluster and a supermassive black hole called Sagittarius A|$^{*}$| (Sgr A|$^{*}$|⁠; e.g. Genzel, Eisenhauer & Gillessen 2010). The combination of these two factors gives rise to a unique, high-energy environment within our Galaxy, which challenges our understanding of key physical processes. The origin of the S-star cluster, for instance, remains unknown to this day, because the strong tidal force in this region inhibits standard star formation from molecular clouds (Genzel et al. 2010). In addition, the initial mass function (IMF) of stars near the GC has been a long-standing point of debate, whose resolution can give us an important clue as to the mass assembly history of the GC. Some studies find it is consistent with the canonical Kroupa (2001) IMF (Maness et al. 2007; Löckmann, Baumgardt & Kroupa 2010), while others find evidence for a top-heavy IMF in disc structures near Sgr A|$^{*}$| (e.g. Paumard et al. 2006; Klessen, Spaans & Jappsen 2007; Bartko et al. 2010; Lu et al. 2013; von Fellenberg et al. 2022). The combination of high line-of-sight extinction and source crowding make the study of the GC challenging, requiring specialized instruments such as GRAVITY (Eisenhauer et al. 2011).

One promising means of better understanding our GC comes from hypervelocity stars (HVSs), which are stars ejected from the GC at extremely high velocities. They are believed to be ejected from the vicinity of Sgr A|$^{*}$| (order of 10 au), but can be observed in parts of the sky more accessible to study (Hills 1988; Yu & Tremaine 2003). This in turn allows for detailed stellar parameter measurements with existing observatories, and the study of wavelength ranges entirely inaccessible near the GC. It has also been suggested that HVSs could explain the existence of the S-star cluster, because the companion of the HVS in the progenitor binary is deposited at radii consistent with the S-star cluster (Gould & Quillen 2003). This migration of already-formed stars into the S-star cluster could therefore solve the issue of in situ star formation in the strong tidal field near Sgr A|$^{*}$| (e.g. Ghez et al. 2003, see, however, Habibi et al. 2017).

In recent years, the term HVS has been used inconsistently. For clarity, we define an HVS to be characterized exclusively by having been ejected from the vicinity of Sgr A|$^{*}$|⁠. Since ejections can occur at a range of velocities, some of these HVSs will still be bound to the Galaxy.

The discovery of HVSs has proved challenging in spite of their significant scientific potential. Only a few dozen promising candidates (e.g. Brown, Geller & Kenyon 2014) and a single star that can be confidently traced back to the GC (Koposov et al. 2020) have been identified. The main challenges include their rarity and the difficulty in disentangling these stars from e.g. hyper runaway stars from the Galactic disc (e.g. Przybilla et al. 2008; Kreuzer, Irrgang & Heber 2020; Irrgang et al. 2021).

In this work, we perform a targeted survey of HVS candidates. The candidates are identified using a novel, highly efficient selection method with Gaia Data Release 3 (DR3; Gaia Collaboration 2016, 2023), which means we have excellent sensitivity to HVSs from relatively few candidates. We perform ground-based follow-up observations for the most promising HVS candidates and use these observations in combination with sophisticated models to constrain the GC environment and the ejection of HVSs. We use these constraints to investigate the total population of HVSs accessible to Gaia in current and future data releases. We furthermore evaluate previously identified candidates using our observational constraints.

In Section 2, we describe how we select HVS candidates from Gaia DR3. In Section 3, we present our follow-up observations. In Section 4, we discuss the model and simulations, in addition to how we apply the observational selections to the models. In Section 5, we present constraints on the GC environment and HVS from our observations in combination with the simulations. In Section 6, we discuss our results and implications for existing HVS candidates as well as what we might expect from Gaia DR4. Finally, in Section 7 we provide closing remarks.

2 HVS CANDIDATE SELECTION

The main challenge when searching for HVSs is that they are extremely rare; there are for instance no HVSs with velocity greater than 700 km s|$^{-1}$| in the 34M sources in Gaia DR3 with radial velocities, that could be robustly identified (Marchetti, Evans & Rossi 2022). HVS candidates selected from observational catalogues for follow-up observations are therefore typically strongly contaminated by non-HVSs. Here, we present our new method to identify HVS candidates originating from the GC. In Section 2.1, we present the general principle we use to find HVSs, and in Section 2.2 we provide the precise selections on the Gaia catalogue for our observational campaign.

2.1 Radial trajectory

Let |$\vec{R}$| be the vector from the GC to the star and |$\vec{V}$| the star’s velocity vector relative to the GC. We can write these down in terms of the observables

• |$\vec{R_0}$|⁠: Vector from the GC to the Sun

• |$\hat{n}$|⁠: Unit vector from the Sun to the source

• D: Distance from the Sun to the source

• |$\vec{V_0}$|⁠: Velocity of the Sun in the Galactic frame of rest

• |$V_{\mathrm{r}}$|⁠: Radial velocity of the source relative to the Sun

• |$\vec{\mu }$|⁠: Proper motion of the source in the Heliocentric frame (always |$\perp$| to |$\hat{n}$|⁠)

giving

and

We provide a sketch of the configuration in Fig. 1. We use the measurement of |$\vec{R_0}$| from GRAVITY Collaboration (2018), which is 8.122 kpc, and |$\vec{V_0}$| from Drimmel & Poggio (2018), which is |$[12.9, 245.6, 7.78]$||$\text{km}\, \text{s}^{-1}$|⁠. The unit vector |$\hat{n}$| is set by the sky coordinates of a star which are provided by Gaia, as is the proper motion vector |$\vec{\mu }$|⁠. For most sources in Gaia DR3, approximately 1.5 billion, we have measurements for all but D and |$V_{\rm r}$|⁠. It is worth noting here that Gaia does measure parallaxes, but only for nearby sources do these result in precise distance measurements. The way of circumventing this lack of knowledge of D and |$V_{\rm r}$| is to only focus on objects that have radial trajectories pointing out from the GC. As we will see later, this enables us to calculate the distance and radial velocity from proper motion and sky position alone. For a star on a radial trajectory from the GC, its velocity aligns with the vector pointing from the GC to the star. This means the position |$\vec{R}$| and velocity |$\vec{V}$| are parallel and thus obey

We can expand this using equations (1) and (2) leaving the terms

and we can pull out the distance and radial velocity by taking the dot products with |$\hat{n}$| and |$\vec{\mu }$| and rearranging.

Taking the dot product of equation (4) with |$\hat{n}$| reveals

and similarly the dot product of equation (4) with |$\vec{\mu }$| gives

This means that for a star moving on a radial trajectory, the distance and radial velocity are set by the observables available for |$\sim 1.5$|B stars in Gaia DR3.

For HVSs, the assumption that their trajectory is purely radial can be reasonable, since stars moving much faster than the escape velocity are not deflected significantly by the Galactic potential (Kenyon et al. 2018; Boubert, Erkal & Gualandris 2020). In this work, we apply the above presented method to determine the distance and radial velocity to all stars in Gaia DR3. Since those solutions are only physical for stars on radial trajectories, which most stars in Gaia are not, we refer to them as the implied distance, |$D_{\rm I}$|⁠, and implied radial velocity, |$V_{\rm r, I}$|⁠.

We convert |$D_{\rm I}$| to implied parallax (⁠|$\varpi _{\rm I}$|⁠) to allow for convenient comparison to measurements by taking the inverse of |$D_{\rm I}$|⁠. We calculate the uncertainty on |$\varpi _{\rm I}$| assuming that the proper motion measurement is the only source of uncertainty. We make the simplifying assumption that the proper motion uncertainties are uncorrelated and since |$\varpi _{\rm I}\propto \mu$|⁠, we can approximate the uncertainty on |$\varpi _{\rm I}$| as

with |$\hat{\mu }{_{\delta }}$| and |$\hat{\mu }_{\alpha }$| the basis vectors for the proper motion in Dec and RA respectively, and |$\sigma _{\mu _\delta }$| and |$\sigma _{\mu _{\alpha }}$| the uncertainties on the proper motions in Dec and RA, respectively. For a more accurate uncertainty analysis one would sample over the uncertainty on all observables.

2.2 HVS candidate selection criteria

The |$\varpi _{\rm I}$| and |$V_{\rm r, I}$| we derived in the previous section can be used to identify HVS candidates. Here, we give a short conceptual overview of our approach. We first compare the implied parallax to the Gaia parallax and reject candidates for which they are inconsistent. In addition, we determine the location in the Hertzprung–Russel (HR) diagram according to the implied distance and remove any candidates for which we consider this solution to be unphysical. Furthermore, we only look at HVS candidates with a high enough 3D implied velocity for our radial trajectory assumption to be reasonable. Lastly, we apply additional cuts that are aimed to reduce the contamination of our sample. We base some of these cuts on catalogues of mock HVSs obtained from simulations described in Section 4.2. These simulations model the orbits of ejected HVS populations and use the characteristics of Gaia to predict their observational properties.

Having discussed the general principle of how we select HVS candidates in Gaia, we now provide the detailed requirements set for a particular star to be observed in our survey. We categorize these into four groups, each discussed in detail below and summarized in Table 1. The goal of the selections presented here is to increase the purity of the sample; i.e. reduce the number of non-HVSs while retaining the real HVSs as much as possible. This is motivated by the necessarily limited number of sources we can provide follow-up observations for.

2.2.1 Astrometric and kinematic selections

We start by using the astrometric and kinematic properties of stars to select potential HVSs. HVSs need to be very fast to prevent significant deviation from their otherwise radial trajectories by the torque on the orbit generated by the non-sphericity of the Galactic potential. For this reason, we only consider stars with total implied velocity, |$V_{\text{I}}$|⁠, in the range |$[800, 3500]$||$\text{km}\, \text{s}^{-1}$|⁠, independent of Galactocentric distance. The lower limit ensures that (excluding the inner GC) all stars would be unbound. The upper limit includes most of the extremely fast predicted HVSs (Rossi, Kobayashi & Sari 2014), while limiting contamination by non-HVSs (which can have extremely high implied velocities). Second, we compare the implied parallax to the parallax measured by Gaia. We only select sources with implied parallax consistent with measured parallax within |$2\sigma$|⁠, considering both the uncertainty on the measured and the implied parallaxes. We do this comparison to the parallax and not the distance since the parallax has Gaussian uncertainties.

The source of uncertainty on |$\varpi _{\text{I}}$| is the measured uncertainty on |$\vec{\mu }$| from Gaia. Because comparing parallaxes is not constraining when both are highly uncertain, we only consider sources where |$\varpi _{\text{I}}$| over its associated uncertainty is larger than five. In addition, most parallaxes in the Gaia catalogue have large fractional uncertainties. This also means that requiring the implied parallax and measured parallax to be consistent is often not a stringent requirement. To alleviate this, we additionally compare the implied distances to the photogeometric ones determined in Bailer-Jones et al. (2021). The photogeometric distances make use of the colour–magnitude information from stars by incorporating a colour–magnitude dependent prior on the extinction-corrected absolute magnitude, which varies as a function of sky position. Similar to the parallax selection, we require that the implied distances are consistent with the photogeometric ones. As a threshold we use |$2\sigma$|⁠, where we use the 16th and 84th percentiles from Bailer-Jones et al. (2021) as the negative and positive |$1\sigma$| estimates for the posterior respectively.

2.2.2 Photometric selections

The implied position of many stars in the HR diagram (according to the implied distance) will be unphysical, because the implied distances are only appropriate for stars moving on radial trajectories. To remove unphysical solutions, we only consider sources within part of the HR diagram.

We first correct the G, |${\rm G_{RP}}$|⁠, and |${\rm G_{BP}}$| magnitudes for extinction using the 2D extinction map from Schlegel, Finkbeiner & Davis (1998) and the recalibration from Schlafly & Finkbeiner (2011), assuming every source is behind the extinction layer. We use an |$R_{\rm V}=A_{\rm V}/{\rm E(B-V)}$| of 3.1 with the Fitzpatrick (1999) extinction law. We use the values provided by Gaia to convert the extinction at 550 nm to the relevant bands.¹ Since those corrections are formally only valid when the intrinsic colour is known, we perform 10 iterations over the extinction coefficients, adjusting the colour of the source at each step. We only select sources which are implied to be within |$\Delta {\rm M_G}=2$| of the main sequence, which we define to be |${\rm M_G} = 4.3\times ({\rm G_{BP}-G_{RP}}) + 0.5$|⁠. This is a fairly simplistic selection, which can be improved upon in future surveys by, for instance, taking the relative contamination of our HVS candidate catalogue by field stars into account as a function of position on the HR diagram. To exclude highly extincted regions, we only consider HVS candidates with a G-band extinction coefficient |$A_{\text{G}}\lt 1.5$|⁠.

We noticed that the relative contamination of our HVS candidate catalogue by field stars was lowest for blue sources, based on the simulations we discuss in Section 4. For this reason we select only sources with an extinction-corrected |${\rm G_{BP}-G_{RP}}$| colour below 0.5, which corresponds to a mass |$\gtrsim 1.3$| M|$_\odot$|⁠. We show these selections in Fig. 2 in combination with a reference HR diagram for sources in Gaia. To provide perspective, we also overlay the density of HVS candidates in Gaia in Fig. A1 in the appendix.

2.2.3 Sky coordinate selections

We know that the distribution of HVSs on the sky should be anisotropic, even for an isotropic ejection mechanism in the GC due to on-sky projection. Moreover, we expect that the contamination of our HVS candidate sample is highest near the GC and anticentre, where our selections are relatively ineffective. In addition, in the Galactic plane we expect the contamination to be higher than at high Galactic latitudes due to the high relative number of field stars to expected HVSs. In order to limit contamination by field stars, we define the expected number of HVSs to number of candidates for different directions in the sky. For this purpose, we use healpix² (⁠|$nside = 3$|⁠; Górski et al. 2005; Zonca et al. 2019) to divide the sky into equal area pixels. We determine the expected density of HVSs across the sky using the simulations described in Section 4. We only observe sources in healpix pixels where the ratio of the density of expected HVSs to candidates (normalized by the number of candidates) is larger than five.³

To prevent contamination from stars in the Large Magellanic Cloud (LMC) or Small Magellanic Cloud (SMC), we additionally require any HVS candidate to be more than 8 and 3 deg from their respective centres.⁴

The selection steps up to this point leave us with only 600 HVS candidates, the previously identified S5-HVS1 being one of the top candidates (Koposov et al. 2020). We provide an overview of the selections in Table 1. The complete source list can be found online.⁵ We provide the first five sources in Table 2.

2.2.4 Instrument-specific selections

We used two observatories for follow-up observations of our candidates: the Isaac Newton Telescope (INT) and the New Technology Telescope (NTT). The INT is a 2.54-m telescope located on the island of La Palma and the NTT is a 3.58-m telescope sited at La Silla in Chile. We used two observatories to gain both Northern and Southern hemisphere coverage of the sky. The instrument-specific selection for the INT is

• |$\textrm {G}\lt 19$|

• |$\textrm {Dec} \gt -20$| deg

• |$\textrm {RA} \lt 60$| deg or |$\textrm {RA} \gt 230$| deg.

Our instrument specific selection for the NTT is

• |$\textrm {G}\lt 19.3$|

• |$\textrm {Dec} \lt 0$| deg.

We find a total of 284 HVS candidates in Gaia that match all of the listed selections. We performed spectroscopic follow-up observations for these sources, which are described in Section 3.

3 OBSERVATIONS

In the previous section, we discussed how we select our HVS candidates. In this section, we describe the follow-up observations we performed for those candidates. We describe the instrument set-up in Section 3.1, the data reduction in Section 3.2, the spectral analysis in Section 3.3, and lastly the observed stars in Section 3.4.

An overview of our observations is given in Table 3.

3.1 Instrument set-up

To confirm or reject HVS candidates, we do not require high radial velocity precision. We therefore set up the instruments to allow efficient radial velocity measurements with uncertainties of a few tens of |$\text{km}\, \text{s}^{-1}$|⁠.

At the INT, we use the Intermediate Dispersion Spectrograph (IDS) with the EEV10 detector and the R400V grating. This gives us an approximate resolution of |$\textrm {R} = \frac{\lambda }{\Delta \lambda } = 1600$| at 4500 Å.⁶ We were allocated a total of 12 nights of observing at the INT.

At the NTT, we use the ESO Faint Object Spectrograph and Camera (v.2; EFOSC2) with grism Gr#07. This set-up gives us a resolution of about 500 at 3800 Å.⁷ We were allocated seven nights of observing at the NTT.

3.2 Data reduction

In the following, we describe the data reduction steps applied to the raw data from the INT and NTT. We use the same pipeline for the two telescopes, with minor adjustments where needed to account for the different instruments.

We perform standard bias subtraction and flat-field correction to all science spectra. We then mask out any cosmic rays using ccdproc (Craig et al. 2022). We trace and subsequently extract the spectra using pyraf, which is the python implementation of iraf (Tody 1986, 1993). Background subtraction is performed by selecting a region around the science spectrum. Through visual inspection, we make sure there are no sources in the background region and make adjustments where needed. In the case of poor seeing during the observation, for instance, the background region might need to be further separated from the centre of the science spectrum. In the case of the INT, wavelength calibration is performed using the CuAr + CuNe lamp. For the NTT, we make use of the He and Ar lamps. An arc spectrum is taken after every science exposure to account for instrument flexures. No flux calibration is applied, since it is not required for radial velocity measurements.

We additionally calculate the signal-to-noise (S/N) as a function of wavelength. We start by taking the square root of the sum of squares of the photon noise and the readout noise. We remark that this is an approximation, because in the regime of low photon counts, the constraint on the flux is not Gaussian (see e.g. Guy et al. 2023, section 4.2.9). We use the same trace used in extracting the science spectrum to extract the noise spectrum and perform error propagation in subtracting the background.

3.3 Spectral analysis

Having obtained wavelength-calibrated spectra, we now describe our analysis. We determine the radial velocities from our reduced data using rvspecfit (Koposov 2019). This software relies on direct pixel fitting using interpolated stellar templates from the PHOENIX spectral library (Husser et al. 2013). For the spectra from the NTT we use the wavelength interval |$3500 \lt \lambda \lt 5240$| Å, while for the INT we use the wavelength interval |$3600 \lt \lambda \lt 6800$| Å. We perform barycentric correction on the resulting radial velocities using astropy.

3.4 Observed HVS candidates

We obtained radial velocity measurements for 196 out of our 284 HVS candidates. This translates to a completeness factor of about 69 per cent. We only consider spectra with an S/N larger than two and three for the INT and NTT, respectively. For lower S/N spectra, we were not able to recover reliable radial velocity measurements. For another five stars, we found previous radial velocity measurements in the SDSS DR14 low-resolution stellar catalogue (Abolfathi et al. 2018) or in the LAMOST DR8 low-resolution catalogue (Zhao et al. 2012). We give the data on these stars in Table B1. We will also consider these five sources in our observational sample, bringing our total completeness to about 71 per cent. In Fig. 3, we plot the sky distribution of our HVS candidate catalogue, indicating which ones have been observed successfully.

We provide an observational log of our observations in an online table,⁸ the first five rows of which are shown in Table 4.

In Fig. 4, we plot the implied radial velocity against the measured radial velocity. We can see that only the previously identified S5-HVS1 falls on the bisection, clearly identifying it as an HVS. However, not all HVSs will be on the bisection, because either their trajectories are not purely radial or due to observational uncertainty. In practice, this means we cannot state with absolute certainty that the other observed sources are not HVSs, because of observational uncertainty and orbits which are not purely radial. We will return to this point in Section 4.3.

4 HVS MOCK CATALOGUES

In this work, we use a comparison between observations and simulations to infer properties of HVSs. Section 4.1 describes the HVS ejection mechanism and Section 4.2 discusses the simulation suite in which it is applied. We use these simulations to make informed decisions for our HVS candidate selection, mentioned in Section 2.2. Most importantly, the simulations allow us to provide constraints on the ejection of HVSs and the GC environment as presented in Section 5.

4.1 Hills mechanism

We assume HVSs are ejected through the Hills mechanism proposed in Hills (1988; for a review see Brown 2015). In this mechanism, a binary star system approaches a massive black hole (MBH) within the tidal radius, where the tidal force from the MBH exceeds the binary self-gravity. The result is that the binary is separated; one star is captured by the MBH in a tight, elliptical orbit and the other star is ejected as an HVS. If the progenitor binary is on a parabolic orbit, which is an excellent approximation in the scenario assumed here (Kobayashi et al. 2012), both stars in the binary have an equal likelihood of being ejected, independent of mass ratio (Sari, Kobayashi & Rossi 2009). The velocity at which the HVS is ejected from the sphere of influence of the MBH is given by

where we neglected a multiplicative constant factor of order unity which depends on the geometry of the encounter (see e.g. fig. 10 in Sari et al. 2009). In equation (8), |$m_{\text{c}}$| is the mass of the captured companion, |$M_{\text{MBH}}$| in the case of the Galaxy is the mass of Sgr A|$^{*}$| (⁠|$4\times 10^6$| M|$_\odot$|⁠; Eisenhauer et al. 2005; Ghez et al. 2008; Gillessen et al. 2009, 2017; Schödel, Merritt & Eckart 2009; Boehle et al. 2016; Do et al. 2019; GRAVITY Collaboration 2019), and M the total mass of the progenitor binary (Sari et al. 2009; Kobayashi et al. 2012; Rossi et al. 2014). Most of the stars ejected through this mechanism are still bound to the Galaxy, but the ejections can occur at velocities in excess of 1000 |$\text{km}\, \text{s}^{-1}$|⁠, putting stars on orbits unbound to the Galaxy.

4.2 Simulations

The simulation framework we use is speedystar⁹ (Contigiani, Rossi & Marchetti 2019; Evans et al. 2022b). This code is able to create mock observational catalogues of Hills mechanism HVSs, and of stars ejected by a binary MBH (Evans et al. 2023). Here, we focus exclusively on the Hills mechanism. Below we describe the most important steps to understand the model, for details we refer to Evans et al. (2022b).

We start by creating an HVS progenitor population of binaries. These binaries are characterized by the primary zero-age main-sequence mass (⁠|$m_{\text{p}}$|⁠), mass ratio (q), and the orbital semimajor axis (a) at the moment of disruption. The mass of the primary is drawn from a mass function (MF) which follows a power-law distribution of the form |$f(m_{\text{p}})\propto m_{\text{p}}^{-\kappa }$| in the range |$[0.1, 100]$| M|$_\odot$|⁠. The secondary is then assigned a mass according to a mass ratio distribution of the form |$f(q)\propto q^\gamma$|⁠, with |$0.1 \le q \le 1$|⁠. Lastly, the orbital semimajor axis are assigned assuming a binary orbital period distribution of the form |$f(\log P) \propto (\log P)^\pi$|⁠. Different values for the three parameters |$\kappa$|⁠, |$\gamma$|⁠, and |$\pi$| are explored with uniform priors to account for uncertainty in the HVS progenitor population. For |$\kappa$| we use the interval |$[0.3, 3.3]$| to allow for both very top-heavy and bottom-heavy MFs. For |$\gamma$| we use the interval |$[-2, 2]$| and for |$\pi$| the interval |$[0, 2]$|⁠. We allow for the wide range in MF index because of the uncertainty on the IMF near the GC, alluded to in the introduction. IMF indexes of 0.45, 1.7, and a canonical 2.3 have all been claimed near Sgr A|$^{*}$| (Bartko et al. 2010; Löckmann et al. 2010; Lu et al. 2013, respectively). The parameter range on |$\pi$| allows for a log-uniform period distribution (also known as Öpik’s law), or distributions with increasing fractions of wide binaries. The parameter range on q encapsulates the results from different studies into massive binaries in star-forming regions (Sana et al. 2012, 2013; Moe & Stefano 2013, 2015; Dunstall et al. 2015 ).

For each star we assume solar metallicity, with |$Z_\odot =0.0142$| (Asplund et al. 2009). Additionally, we performed tests with supersolar metallicity values of |$[\mathrm{Fe/H}]=0.3$|⁠, since stars near the GC are mostly metal rich (Do et al. 2015; Feldmeier-Krause et al. 2017; Nandakumar et al. 2018; Schultheis et al. 2019).

From this progenitor population, binaries are drawn at a constant rate. One of the binary components is selected at random and assigned an ejection speed according to equation (8). We assume the ejection is isotropic and assign the star a random age within the total lifetime of the shorter lived star of the progenitor binary, where we implicitly assume a constant star formation rate. The ejected stars are subsequently integrated in the potential of the Galaxy (see Marchetti et al. 2018, table 1). The result is a simulated current population of HVSs throughout the Galaxy. These HVSs are evolved using single stellar evolution models from Hurley, Pols & Tout (2000) in AMUSE (Portegies Zwart et al. 2009; Pelupessy et al. 2013; Portegies Zwart et al. 2013; Portegies Zwart & McMillan 2018).

Lastly, mock observations are performed incorporating the visual extinction along the line of sight, effective temperature, surface gravity, and metallicity of the simulated HVSs. This is done using the MESA Isochrone and Stellar Tracks model grid (Choi et al. 2016; Dotter 2016). This allows us to determine the magnitude and colour of these stars as they would be observed by various observatories. Astrometric uncertainties in parallax and proper motion, as observed by Gaia, are determined using the DR3 astrometric spread function of Everall et al. (2021).

4.3 HVS candidate selection applied to simulations

Now that we have described both our simulations and observations, we will describe how we apply the selections discussed in Section 2 to our simulated HVSs.

First, we have to consider which stars are observed by Gaia in the first place. Although the magnitude limit of Gaia is |$\text{G}=20.7$| (Gaia Collaboration 2016), not all sources down to that depth appear in the published Gaia catalogue for various reasons. A description of which stars do and do not appear in the published Gaia catalogue is referred to as the Gaia selection function, which has been empirically determined in Cantat-Gaudin et al. (2023).¹⁰ Although Gaia is virtually complete down to our magnitude limit of |${\rm G}=19.3$| except in crowded regions, we none the less incorporate the empirical Gaia selection function for completeness. We also convolve the ‘true’ parallaxes and proper motions with their uncertainties. No uncertainties are considered on the photometric measurements.

Most of the selections described in Section 2.2 can be directly applied to the resulting mock HVS catalogue. We can calculate the implied distance and radial velocity for the mock HVSs in the same way we did for the Gaia catalogue and apply our selections. The same can be done for the photometric, sky coordinate, and instrument-specific selections. An exception is the photogeometric distance from Bailer-Jones et al. (2021). To apply this selection to the mock HVS catalogue, we resort to the Bailer-Jones (2015) source code (C.A.L. Bailer-Jones, private communication) and calculate the photogeometric distances for the simulated HVSs. After obtaining the photogeometric distances for all our mock HVSs in this way, we again apply the same selections to our mock catalogue as to the Gaia catalogue. The only exception is the RUWE, which we do not model for our simulated HVSs. We expect HVSs to be single stars and for single stars fewer than one in a million have |${\tt RUWE}\gt 1.4$| (Penoyre, Belokurov & Evans 2022).

Since our observational catalogue does not cover every candidate from the selections introduced in Section 2, we need a prescription to quantify our observational completeness. In Fig. 5, we show our follow-up observational completeness as a function of the apparent magnitude of the HVS candidates. For each bin in G magnitude, we calculate the fraction of sources we observed over the number of sources in our catalogue (see Section 2.2.4). As can be seen, the observational completeness decreases for fainter sources, since we prioritized bright sources. In addition, we have lower completeness for regions with many HVS candidates that were not accessible for a significant amount of time during the observational runs. This expresses itself mainly in a lower completeness near the GC. In Fig. 6, we show our follow-up observational completeness as a function of the angle to the GC. Since the observational completeness shows a stronger bias in G magnitude than in the sky position, we decide to model the observational completeness only on G magnitude. We bin the magnitudes of our candidates using 50 bins between |${\rm G}=0$| and |${\rm G}=19.4$| and calculate the fraction of observed stars within each bin. If there are no candidates in a bin, we consider the observations to be complete. The model of the selection function of our observations is then the completeness fraction per bin in G magnitude, which we can apply to our simulations.

Lastly, we noted in Section 3.4 that we cannot reject with absolute certainty the HVS candidates that have significantly different observed and implied radial velocities. To account for this, we only consider simulated HVSs for which |$V_{{\rm r,I}}-V_{{\rm r}} \lt 300$||$\text{km}\, \text{s}^{-1}$|⁠, which holds for about 86 per cent of HVSs in our simulations. None of the sources we observed, except S5-HVS1, meet this criterion.

By combining all observational selections and applying them to the simulations, we can predict the number of HVSs that should be in our observational sample, which we use to provide constraints on the GC in Section 5.2.

5 CONSTRAINTS ON THE GALACTIC CENTRE

Having described the observations, we can now provide physical constraints on the GC environment and HVS ejection process. Before we do so in Section 5.2, we will first discuss the properties of HVSs that can be observed with Gaia in Section 5.1.

5.1 HVSs observable by Gaia

Not all HVSs can be observed by Gaia. Low-mass stars, for instance, can only be observed out to short heliocentric distances, due to their intrinsic low brightness. To establish the properties of HVSs that can be observed by Gaia, we make use of the simulations described in Section 4.2. We require that to be observed by Gaia and be considered an HVS, the star must have an apparent G magnitude below 20.7, in addition to travelling on an unbound orbit. We additionally consider the Gaia selection function mentioned in Section 4.3.

For every set of model parameters we can now find the stars that would appear in the Gaia catalogue. We caution that due to observational uncertainties, not all of these real HVSs could be trivially identified as such; in practice, it can be challenging to determine if a star is in fact unbound or not. In Fig. 7, we plot the distance against the mass for HVSs that can be observed by Gaia, along with the colour–magnitude diagram. In the left column, we add the contributions from all MFs within our prior range (see Section 4.2). In the middle and right columns, we show the distribution for two specific MFs (a top-heavy one, |$\kappa =1.7$|⁠, and one close to the canonical Salpeter MF, |$\kappa =2.3$|⁠). The resulting distribution should be seen as a prediction for the overall population distribution of HVSs present in the Gaia catalogue, discovered or undiscovered. We can see that most of the HVSs in the Gaia catalogue will be at heliocentric distances between about 10 and 200 kpc, be of the order of a few solar masses, and be main-sequence stars. Below heliocentric distances of |$\sim 5$| kpc, the volume covered by Gaia is too small which makes discovering any HVSs within this region highly unlikely. The vertical feature in the top row at |$\sim 8$| kpc is the GC, where, due to projection effects, an overdensity of HVSs on the sky is present. The cut-off feature towards low-mass/distant stars is caused by the detection limit of Gaia in apparent magnitude. The stars below this boundary are giants, which are visible out to further distances than main-sequence stars. The upper right of the figure is not populated, because extremely massive stars do not live long enough to travel far from the GC. In Fig. C1, we additionally show the apparent magnitude distribution for an MF index |$\kappa =2.3$|⁠.

Interestingly, the MF has little effect on the ratio of evolved to main-sequence stars we expect to be within the Gaia catalogue. The fraction of evolved HVSs is expected to be about 0.17, changing by less than 0.03 between different MFs.

Using these results, we can also determine the completeness of Gaia for HVSs as a function of mass. In other words, what fraction of HVSs is expected to be in the Gaia catalogue as a function of mass. For this, we only consider HVSs that have not evolved into remnants. We show the results in Fig. 8 for varying distance limits. We can see that the completeness for HVSs in Gaia approaches zero for masses below |$\sim 1 \,\, {\rm M_\odot }$| even when only considering HVSs within 50 kpc from the Sun, which are the ones our observations are mainly sensitive to. For high-mass HVSs (M|$\gtrsim 5$| M|$_{\odot }$|⁠), the completeness converges to nearly one, because those stars do not live long enough to travel far from the GC.

5.2 HVS parameter constraints

Having established the types of HVSs that Gaia is sensitive to, we will now describe our constraints on the ejection of HVSs.

By applying the same selections to the simulations described in Section 4 as the observations, we can predict the number of HVSs that should be observed in our catalogue given a set of model parameters. By rejecting models that are incompatible with the observed number of HVSs, we provide constraints on the model parameters, and as such, the physical environment of the GC and the process of HVS ejection.

We use uniform priors within the parameter ranges given in Section 4.2. The posterior probability of the model with parameters M is therefore given by

with |$N_{\rm HVS}$| the number of observed HVSs and |$\lambda$| the average number of expected HVSs according to the simulations for a given set of parameters M.

We find that the main parameters we can constrain based on the number of observed HVSs are the ejection rate and MF (as is also the case in Evans et al. 2022a, b). Both the log-period distribution and mass-ratio distributions are significant in their effect on the observed number of HVSs, but affect the results much less than the MF and ejection rate.

In Fig. 9, we provide the posterior on the MF and ejection rate, marginalized over the binary period and mass ratio distributions. Compared to previous studies, we see a remarkable improvement in the upper limit on the combined ejection rate and MF. Even though we only obtained 196 radial velocity measurements, our results are about 1 dex more constraining than the |$\sim 34$| M radial velocity measurements in Gaia DR3 (also considering S5-HVS1 was discovered). This demonstrates the effectiveness of our targeted survey versus blind surveys when identifying HVSs. It is also worth noting that we use a more conservative prior on the log-period distribution compared to earlier work, since the prior used in Evans et al. (2022a, b) results in higher ejection velocities to which our observations are most sensitive. We only highlight the upper limit and not the lower limit, because the upper limit is better constrained. The lower limit is highly sensitive to the prior used (e.g. uniform or Gamma distribution). In Fig. D1, we show the posterior if we use a log-uniform prior of |$1/\lambda$|⁠. The posterior with the log-uniform prior shows that the upper limit remains, while allowing lower ejection rates with equal probability.

For |$\kappa \gtrsim 1.5$|⁠, we can see in Fig. 9 that the ejection rate is a strong function of the MF index. This is caused by our limited sensitivity to low-mass stars, as demonstrated in Fig. 7. A more bottom-heavy MF would cause a lower fraction of HVSs to be observable by Gaia, thus necessitating a relatively higher ejection rate. For extremely top-heavy MFs, the ejection rate and MF index become relatively independent due to two competing effects:

• a larger fraction of HVSs is sufficiently massive to be observed with increasingly top-heavy MFs, and

• the finite lifetime of these increasingly massive stars causes a significant fraction of HVSs to evolve into stellar remnants before exiting the sphere within which they would be visible were they main-sequence stars.

Overall, our upper limit ejection rate constraints are consistent with the process where stellar binaries’ centre of mass trajectories have undergone gravitational scatterings off other (single) stars or giant molecular clouds: once every |$10^{4}-10^{6}$| yr a binary star is sent on a highly eccentric orbit such that its pericentre becomes equal or smaller than the tidal separation radius (a.k.a. ‘loss cone orbit’) and an HVS is ejected (Yu & Tremaine 2003; Perets, Hopman & Alexander 2007). For a recent review on rates see Stone et al. (2020, especially section 3.6.3). Robust rate predictions to compare with our observational constraints should include realistic descriptions of the phase space and properties of binary stars, as well as those of single stars and giant molecular clouds out to |$\sim 100$| pc from SgrA* in our GC. Such calculations and comparisons are beyond the scope of this paper and will be presented in a follow-up investigation.

We can re-parametrize the ejection rate, such that it is no longer a function of the MF index by only considering HVSs massive enough to be observable by Gaia. The total ejection rate defines the ejection rate of HVSs in the mass range |$[0.1, 100]\ \mathrm{M}_\odot$|⁠, as explained in Section 4.2. We can transform this to the ejection rate in a different mass range by integrating the MF, which allows us to determine the number of stars in a mass interval |$[M_1, M_2]$|

The ejection rate for stars in the mass range |$[M_1, M_2]$| is then the total ejection rate times the fraction of stars in the mass interval |$[M_1, M_2]$|⁠. In Fig. 10, we give the resulting posterior for HVSs more massive than |$1\mathrm{M}_\odot$|⁠, making them potentially accessible to Gaia. The figure demonstrates that irrespective of the MF, the ejection rate of HVSs more massive than |$1\mathrm{M}_\odot$| cannot be higher than about |$10^{-5}$| yr|$^{-1}$| at |$2\sigma$| confidence.

In Fig. 11, we show the flight time (i.e. time since disruption) distribution of HVSs that our observations are sensitive to. We can see that the observations presented in this work are sensitive to ejections over the past |$\sim 50-100$| Myr. Our ejection rate constraint therefore applies to ejections within the last |$\sim 50-100$| Myr. HVS ejections that occurred more than 100 Myr ago would not produce any observable HVSs in our sample and can therefore not be constrained with our observations.

In the hypothetical case we observed two HVSs instead of one, the effect on our constraints is mainly that the posterior (Figs 9 and 10) becomes more narrow (the |$2\sigma$| confidence interval by |$\sim 0.5$| dex). In addition, the |$2\sigma$| upper limit constraint would increase by about |$\sim 0.1$| dex.

5.3 How many HVSs are in Gaia?

Using our constraints from Section 5.2, we can predict how many HVSs are in the complete Gaia DR3 catalogue. Along with the predicted mass–distance relation, shown in Fig. 7, this gives us a perspective on the discovery space of HVSs in Gaia. In Fig. 12, we show the expected number of HVSs in Gaia as a function of the MF. Within |$1\sigma$|⁠, we expect there to be between about 5 and 45 HVSs (⁠|$1.7 \lt \kappa \lt 2.3$|⁠), with the mode of the distribution around 18. The significant uncertainty in this number is caused by the single confirmed HVS within our observational sample. Importantly, the number of predicted HVSs is about one to two orders of magnitude less than predicted in earlier work (Marchetti et al. 2018). This difference is mostly caused by a much lower ejection rate than previously expected. Theoretical prospects relying on Gaia discovered HVSs will have to take this into account.

6 DISCUSSION

Having presented our main results, we will now provide some discussion. We start by reviewing our assumptions in Section 6.1, second we discuss the influence of the LMC on our results in Section 6.2, then we discuss the previously discovered HVS S5-HVS1 in light of our results in Section 6.3, we discuss some previously identified HVS candidates in Section 6.4, and lastly provide prospects for the discovery of additional HVSs in Section 6.5.

6.1 Assumptions

In this work, we rely on simulations to interpret our observational results. Evaluating the assumptions the simulations rely on is therefore important. Most of these have been discussed in Evans et al. (2022b), including alternative MF and mass ratio parametrizations. We focus here on differences.

The MF we use in the simulations is different from the IMF for HVS progenitor binaries. The simulation assigns a random ejection time within the lifetime of a star, which makes the MF we constrain in this work more analogous to the present-day MF. In addition, the use of a single, time-independent MF in the simulations implicitly assumes a constant star formation rate. This is a simplifying assumption, which we know does not hold in detail (e.g. Nogueras-Lara et al. 2020).

Regarding the log-period distribution, we use a more stringent prior compared to Evans et al. (2022b). Power-law indexes in |$\alpha$| (see Section 4.2) below zero would result in higher ejection velocities because of the lower binary separation (see equation 8). Since our observations are most sensitive to the fastest HVSs, our prior is more conservative.

We assumed a solar metallicity for the progenitor binaries throughout this work. We evaluate the effect of metallicity on our results by running identical simulations with |$[{\rm Fe/H}]=0.3$|⁠. The difference on the posterior presented in Fig. 9 is about 7 per cent of the statistical uncertainty. We therefore conclude that metallicity is currently not an important consideration for our work.

We also tested the influence of the adopted Galactic potential on our results. To compare, we ran the simulations for the potential determined in McMillan (2017). This potential is a more massive commonly used potential compared to our default one. The difference on our posterior from Fig. 9 changes by |$0.01-0.1$| dex, which is well below the statistical uncertainty.

6.2 Influence of the Large Magellanic Cloud

Our search for HVS candidates is based on pure radial trajectories. In reality, the non-spherical potential of the Galaxy deviates HVSs from exclusively radial trajectories. Because we compare our observations with simulations that are performed with a realistic, non-spherical potential and an identical selection procedure, these deviations from radial trajectories are not a concern for our results. However, the LMC is not incorporated in the potential model of the Galaxy in which the simulated HVSs are propagated. This might be a concern if the LMC significantly deviates HVSs from their trajectories (see Kenyon et al. 2018; Boubert et al. 2020). In that case, simulations overpredict the number of HVSs we would be able to find and therefore bias the ejection rate constraint towards values which are too small.

We can evaluate the influence of the LMC by a comparison between integrating orbits of simulated HVSs with and without the LMC. We reran the simulations for |$\kappa =2.3$| including the LMC with a mass of |$1.5\times 10^{11}$| M|$_\odot$|⁠, which is near the total mass most studies find and close to half the mass of the Galaxy within 50 kpc (Erkal et al. 2019; Vasiliev, Belokurov & Erkal 2020; Shipp et al. 2021; Koposov et al. 2023). In the run with the LMC, we found an increase of about 2.8 per cent in the number of recovered HVSs in the simulation, which is about half the Poisson noise for those simulations. Our simulation with the LMC did not include the reflex motion of the Galaxy, but since this effect is of the same order as the deflection (Boubert et al. 2020), we conclude that the LMC has no discernible influence on our results.

6.3 S5-HVS1

The serendipitous discovery of S5-HVS1 in the S|$^5$| survey (Li et al. 2019) appears curious. It is the only HVS that can be unambiguously traced back to the GC and is, by far, the fastest out of the stars that are suggested to originate from the Hills mechanism, in addition to being the closest (Brown et al. 2014, 2018; Koposov et al. 2020). An interesting question to pose is then: What is the chance that the S|$^5$| survey included an HVS? To understand this likelihood, we again make use of the simulations described in Section 4.2.

We consider the footprint of the S|$^5$| survey to be any stars with Gaia parallax |$\varpi \lt 3\sigma _\varpi + 0.2$|⁠, mock DECam (Dark Energy Survey Collaboration 2016) photometry |$15 \lt g \lt 19.5$|⁠, |$-0.4 \lt (g-r) \lt 0.1$|⁠, and falling within 2 deg from an S|$^5$| pointing (see Li et al. 2019, table 2) following Evans et al. (2022b). Only HVSs with heliocentric radial velocities of |$V_{\text{r}} \gt 800$||$\text{km}\, \text{s}^{-1}$| are considered, since these were inspected in Koposov et al. (2020). We find that within the |$1\sigma$| posterior probability on the MF and ejection rate from Fig. 9, the number of expected HVSs in S|$^5$| is about |$0.14_{-0.11}^{+0.36}$|⁠. Although unlikely, the discovery of S5-HVS1 in S|$^5$| is consistent with our constraints on the ejection rate of HVSs.

The radial velocity we measure for S5-HVS1 is |$995\pm 12$||$\text{km}\, \text{s}^{-1}$| when considering the wavelength range 3850–5240 Å. We used this wavelength range in particular because the spectrum we obtained for S5-HVS1 was much higher S/N than the other sources, and we noticed the wavelength calibration below 3850 Å showed systematics. This is most likely caused by the lack of calibration lines being available for that part of the spectrum. Our measurement is consistent with the earlier found |$1017\pm 2.7$||$\text{km}\, \text{s}^{-1}$|⁠, so we find no indication of binarity.

Throughout this work, we assume that S5-HVS1 has been ejected by the Hills mechanism. Other mechanisms have, however, not been ruled out. A three-body interaction between a hypothetical intermediate mass black hole and Sgr A|$^{*}$|⁠, for instance, could have ejected S5-HVS1 and is consistent with the S-star eccentricity distribution (Generozov & Madigan 2020). However, Evans et al. (2023) find that the lack of further detected HVSs in Gaia would necessitate a very specific configuration of any hypothetical intermediate mass black hole and Sgr A|$^{*}$|⁠.

6.4 Distant HVS candidates

In general, the predicted HVSs in Gaia (see Fig. 7) appear similar to those found in Brown et al. (2014), being dominated by distant stars of a few M|$_\odot$|⁠. The difficulty with these stars is that due to their distance, their orbits cannot be confidently traced back to the GC (Irrgang, Kreuzer & Heber 2018; Kreuzer et al. 2020). Confusion with other types of sources, such as disc runaways and halo stars is therefore an issue (e.g. Brown et al. 2018).

We use the posterior on the ejection rate from Fig. 9 to make a prediction on the number of HVSs within the MMT HVS Survey (Brown et al. 2014). We reproduce the same colour selections outlined in section 2.1 of Brown et al. (2012). For simplicity, we consider the footprint of SDSS to be |${\rm Dec}\gt 0$| deg for the Northern Galactic cap, and |${\rm Dec} \gt -15$| deg for the Southern Galactic cap with a Galactic latitude |$|l|\gt 30$| deg (for the sky coverage see Aihara et al. 2011). Brown et al. (2018) only study stars with a heliocentric radial velocity transformed to the Galactic frame greater than 275 |$\text{km}\, \text{s}^{-1}$|⁠, which we follow. In their study, Brown et al. (2018) identify seven stars with a probable origin in the GC. We use the posterior from Fig. 9, in combination with the above selections to calculate how many HVSs there are predicted to be within this footprint. We sample the MF index in the range |$[1.7, 2.3]$| and multiply each sample by the posterior on the ejection rate for that MF index. We find that within a 68 per cent confidence interval, we expect there to be between 0.3 and 2.8 HVSs in the MMT survey footprint. For the 95 per cent confidence interval, we expect there to be between 0.1 and 5.1 HVSs respectively. It is thus likely that a number of the stars identified in Brown et al. (2018) are genuine HVSs. Interestingly however, there is thus a |$\gt 2\sigma$| tension between our results and the number of GC origin HVSs found in Brown et al. (2018). Apart from a statistical deviation, we have two main hypotheses as to the origin of this tension: the ejection rate of HVSs has not been constant over the flight time of the stars identified in Brown et al. (2018), or (a fraction of) the stars identified in Brown et al. (2018) have a different origin (supported by Generozov & Perets 2022). To investigate the first hypothesis, we show the distribution of flight times for predicted HVSs in the MMT footprint in Fig. 13. By comparison to Fig. 11, we can clearly see that the simulated HVSs matching the MMT survey selections tend to have been ejected longer ago relative to those our survey presented in this work is sensitive to. The unbound HVS candidates reported to have a probable origin in the GC have flight times between about 50 and 150 Myr (see Brown et al. 2014, table 1), which is consistent with our predictions.

Brown et al. (2018) estimated the ejection rate of unbound HVSs in the mass range |$2.5 \le M_\odot \le 4$| to be |$1.5\times 10^{-6}$| yr|$^{-1}$|⁠. Using our observations, we can provide an independent constraint on the ejection rate in this mass range. We do this following the same procedure as for Fig. 10, but now with the above-mentioned mass range. Our constraints on this ejection rate are relatively insensitive to the MF varying by a factor of about 3. The |$2\sigma$| upper limit peaks at about |$1.3\times 10^{-6}$| yr|$^{-1}$| for an MF index |$\kappa \sim 2.5$|⁠, while the mode of the posterior peaks at about |$2.8\times 10^{-7}$| yr|$^{-1}$| (at the same MF index). Our maximum likelihood prediction is thus more than five times lower compared to the findings in Brown et al. (2018).

6.5 Prospects

Gaia DR4 will include radial velocities for sources with a limiting magnitude of G|$_{\mathrm{RVS}}\lt 16.2$| (Katz et al. 2023). If we consider all unbound HVSs with our ejection rate constraints from Section 5, we predict there to be about |$3^{+4}_{-2}$| HVSs within the Gaia DR4 radial velocity catalogue if the MF index is between 1.7 and 2.3. It is, however, unlikely that all hypothetical HVSs in the footprint of the radial velocity catalogue of Gaia DR4 could be recognized as such. If we additionally require an accurate parallax measurement (⁠|$\varpi /\sigma _{\varpi }\gt 5$|⁠), we only expect there to be |$1^{+1}_{-1}$| HVSs. If on top of that, we only consider HVSs with radial velocities |$\gt 500$||$\text{km}\, \text{s}^{-1}$|⁠, we are left with fewer than one expected HVS.

A potentially promising means to extend to fainter stars is using radial velocities from the low-resolution Gaia|${\rm G_{BP}/G_{RP}}$| spectra (Verberne et al. 2024). These spectra will be published for all Gaia sources in DR4. The radial velocity measurements from these spectra are most accurate for blue sources (⁠|${\rm G_{BP}-G_{RP} \lesssim 0.7}$|⁠), which is the expected main discovery space for HVSs in Gaia (see Fig. 7). Because the nominal uncertainties are high for radial velocity measurements using the |${\rm G_{BP}/G_{RP}}$| spectra, we only consider HVSs with radial velocities greater than 1000 |$\text{km}\, \text{s}^{-1}$|⁠, in addition to only selecting stars with |${\rm G_{BP}-G_{RP}}\lt 0.7$|⁠. Within this parameter range, we expect there to be |$7^{+9}_{-5}$| HVSs in Gaia DR4. However, significant improvements would be required in the technique of obtaining these radial velocities from the |${\rm G_{BP}/G_{RP}}$| spectra to reliably identify these HVSs (Verberne et al. 2024).

In addition to Gaia DR4, large upcoming spectroscopic survey instruments such as DESI (Cooper et al. 2023), WEAVE (Dalton et al. 2014), and 4MOST (de Jong et al. 2019) provide the opportunity to observe much larger numbers of sources. As part of the low-resolution high-latitude survey, WEAVE will observe about |$20\, 000$| HVS candidates in the Northern hemisphere. This is two orders of magnitude larger than the observational sample described in this work. Moreover, WEAVE will allow us to observe stars down to the detection limit of Gaia. The precise survey strategy for WEAVE will be discussed in upcoming work.

7 CONCLUSION

In this work, we presented a novel selection method of HVS candidates. Utilizing ground-based follow-up observations, we were able to reject all newly identified candidates. By combining this null-detection of new HVSs with sophisticated simulations, we were able to significantly improve the constraints on both the ejection rate of HVSs and the MF of the HVS progenitors. We show that the MF and ejection rate are degenerate because Gaia is only sensitive to HVSs more massive than about |$1\mathrm{M_\odot }$|⁠. We are therefore able to provide robust constraints on the ejection rate for HVSs more massive than |$1\mathrm{M_\odot }$| independent on the MF. Over the past |$\sim 50-100$| Myr, the average ejection rate for these stars cannot have been higher than about |$10^{-5}$| yr|$^{-1}$| at |$2\sigma$| significance.

We use our constraints on the ejection rate and MF to evaluate previously identified HVS candidates. We find evidence that either the ejection rate of HVSs has varied significantly over the past |$\sim 150$| Myr, or these previously identified candidates include stars that do not originate in the GC.

Our constraints predict that within a 68 per cent confidence interval, there are between about 5 and 45 unbound HVSs in the complete Gaia catalogue. The majority of these stars will be main-sequence stars with a mass of a few |${\rm M_\odot }$| and be at heliocentric distances between 10 and 100 kpc.

ACKNOWLEDGEMENTS

The authors would like to thank C.A.L. Bailer Jones for kindly making the source code to his work available and Warren R. Brown for useful comments and suggestions. This work is based on service observations made with the INT (programme ING.NL.22B.002) operated on the island of La Palma by the Isaac Newton Group of Telescopes in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. Based on observations collected at the European Southern Observatory under ESO programme(s) 110.23SU and 111.24MP. EMR acknowledges support from European Research Council (ERC) grant number: 101002511/project acronym: VEGA_P. SK acknowledges support from the Science & Technology Facilities Council (STFC) grant ST/Y001001/1. TM acknowledges a European Southern Observatory (ESO) fellowship. FAE acknowledges support from the University of Toronto Arts & Science Postdoctoral Fellowship program and the Dunlap Institute. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This paper made use of the Whole Sky Database (wsdb) created and maintained by Sergey Koposov at the Institute of Astronomy, Cambridge with financial support from the Science & Technology Facilities Council (STFC) and the European Research Council (ERC). This research or product makes use of public auxiliary data provided by ESA/Gaia/DPAC/CU5 and prepared by Carine Babusiaux. For the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.

Software: numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013, 2018, 2022), dustmaps (Green 2018), healpy (Górski et al. 2005; Zonca et al. 2019), sqlutilpy.¹¹ This research made use of ccdproc, an astropy package for image reduction (Craig et al. 2022).

DATA AVAILABILITY

The Gaia data underlying this article were accessed from https://gea.esac.esa.int/archive/. The HVS candidate catalogue used in this work and the results from our analysis of the calibrated spectra can be accessed from https://zenodo.org/doi/10.5281/zenodo.12179452.

Footnotes

REFERENCES

APPENDIX A: HR DENSITY OF HVS CANDIDATES

As mentioned in Section 2.2.2, the implied position for a star that is not moving on a radial trajectory can place it in a unphysical position in the HR diagram. In Fig. A1, we show where HVS candidates in Gaia that match all but our colour–magnitude selections end up in the HR diagram.

APPENDIX B: RADIAL VELOCITIES FROM SDSS AND LAMOST

In Table B1, we list the five sources for which we found literature radial velocity measurements in either SDSS or LAMOST.

APPENDIX C: G MAGNITUDE DISTRIBUTION Gaia OBSERVABLE HVSS

In Fig. C1, we show the apparent G magnitude distribution for predicted HVSs that are observable by Gaia assuming an MF-index|$\kappa =2.3$|⁠. The two vertical features correspond to the main-sequence and red giant branch. We can see that most HVSs in Gaia will be faint, as expected, since the observed volume increases towards faint magnitudes.

APPENDIX D: LOWER LIMIT ON THE EJECTION RATE AND MF

The lower limit on the ejection rate and MF is sensitive to the prior used. Throughout the main text, we used a uniform prior in |$\lambda$|⁠. In Fig. D1, we instead show the posterior if we use a log-uniform prior equal to |$1/\lambda$|⁠. This prior might be more appropriate, since we do not know the magnitude of the rate, but do know it should be positive. This prior effectively removes the lower limit, since the posterior peaks for the limit of the ejection rate approaching zero.