# Hypervelocity star observations constrain the Galactic Centre

Sill Verberne,1 Elena Maria Rossi,1 Sergey E. Koposov,2,3,4 Tommaso
Marchetti,5 Konrad Kuijken,1 Zephyr Penoyre,1 Fraser A. Evans,6,7 Dimitris
Souropanis,8,9 and Clár-Bríd Tohill,8,10
1Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, the
Netherlands
2Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
3Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
4Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK
5European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching
bei Munchen, Germany
6David A. Dunlap Department of Astronomy and Astrophysics, University of
Toronto, 50 St. George Street, Toronto, ON M5S 3H4, Canada
7Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto, ON M5S 3H4, Canada
8Isaac Newton Group of Telescopes, Apartado 321, E-38700 Santa Cruz de La
Palma, Canary Islands, Spain
9Institute of Astrophysics FORTH, 71110 Heraklion, Greece
10University of Nottingham, School of Physics & Astronomy, Nottingham, NG7
2RD, UK E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### 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 use 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⊙ 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 many previously identified HVS candidates.

###### keywords:

Galaxy: centre – Galaxy: nucleus – stars: kinematics and dynamics – binaries:
general – surveys

††pubyear: 2024††pagerange: Hypervelocity star observations constrain the
Galactic Centre–D

##  1 Introduction

The Galactic Centre (GC) is a highly complex environment, hosting both a
Nuclear Star Cluster and a super-massive black hole called Sagittarius A∗ (Sgr
A∗; e.g. Genzel et al., 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 hisotry of the GC. Some studies find it
is consistent with the canonical IMF from Kroupa (2001) (Maness et al., 2007;
Löckmann et al., 2010), while others find evidence for a top-heavy IMF in disc
structures near Sgr A∗ (e.g. Paumard et al., 2006; Klessen et al., 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 specialised 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 characterised 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 et
al., 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 et al.,
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 et al., 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 et al., 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

$\vec{R}=\vec{R_{0}}+D\hat{n}$ (1)

and

$\vec{V}=\vec{V_{0}}+V_{\rm r}\hat{n}+D\vec{\mu}.$ (2)

We provide a sketch of the configuration in Fig. 1.

Figure 1: A diagram of the position ($\vec{R}$) and velocity ($\vec{V}$), in
blue, of an HVS relative to the GC. This can be expressed in terms of a
combination of the star’s observed position and proper motion (pink) and
measurements of the Sun’s position and velocity (black) along with the
distance to the star ($D$) and its radial velocity ($V_{r}$).

We use the measurement of $\vec{R_{0}}$ from GRAVITY Collaboration et al.
(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

$\vec{R}\times\vec{V}=\vec{0}.$ (3)

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

$\vec{R_{0}}\times\vec{V_{0}}+V_{{\rm
r}}(\vec{R_{0}}\times\hat{n})+D(\vec{R_{0}}\times\vec{\mu}+\hat{n}\times\vec{V_{0}})+D^{2}(\hat{n}\times\vec{\mu})=\vec{0}$
(4)

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

$D=\frac{\hat{n}\cdot\left(\vec{R_{0}}\times\vec{V_{0}}\right)}{\vec{R_{0}}\cdot\left(\hat{n}\times\vec{\mu}\right)}=\frac{V_{0}}{\mu}\left(\frac{\hat{n}\cdot\left(\hat{R}_{0}\times\hat{V}_{0}\right)}{\hat{R_{0}}\cdot\left(\hat{n}\times\hat{\mu}\right)}\right)$
(5)

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

$\displaystyle V_{{\rm r}}$
$\displaystyle=\frac{D\vec{V_{0}}\cdot\left(\hat{n}\times\vec{\mu}\right)-\vec{\mu}\cdot\left(\vec{R_{0}}\times\vec{V_{0}}\right)}{\vec{R_{0}}\cdot\left(\hat{n}\times\vec{\mu}\right)}$
(6)
$\displaystyle=V_{0}\left(\frac{\frac{D}{R_{0}}\hat{V}_{0}\cdot\left(\hat{n}\times\hat{\mu}\right)-\hat{\mu}\cdot\left(\hat{R}_{0}\times\hat{V}_{0}\right)}{\hat{R}_{0}\cdot\left(\hat{n}\times\hat{\mu}\right)}\right).$

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 et al., 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\propto\mu$, we can approximate the uncertainty on $\varpi$ as

$\sigma_{\varpi}=\frac{\sqrt{\left[\vec{R}_{0}\cdot\left(\hat{n}\times\hat{\mu}_{\delta}\right)\sigma_{\mu_{\delta}}\right]^{2}+\left[\vec{R_{0}}\cdot\left(\hat{n}\times\hat{\mu}_{\alpha*}\right)\sigma_{\mu_{\alpha*}}\right]^{2}}}{\hat{n}\cdot\left(\vec{R}_{0}\times\hat{V}_{0}\right)},$
(7)

with $\hat{\mu}{{}_{\delta}}$ and $\hat{\mu}_{\alpha*}$ the basis vectors for
the proper motion in declination (Dec) and right ascension (RA) respectively,
and $\sigma_{\mu_{\delta}}$ and $\sigma_{\mu_{\alpha*}}$ the uncertainties on
the proper motions in Dec and RA respectively (where we take
$\sigma_{\mu_{\alpha}*}=\sigma_{\mu_{\alpha}}\cos^{2}\delta$). 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 categorise these into four groups, each discussed
in detail below and summarised 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}$. Secondly, 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 photo-
geometric ones determined in Bailer-Jones et al. (2021). The photo-geometric
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 photo-geometric 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 et al. (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
bands111https://www.cosmos.esa.int/web/gaia/edr3-extinction-law. 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}}<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⊙. We show these selections in Fig. 2 in combination with a
reference HR diagram for sources in Gaia.

Figure 2: Reference HR diagram created from sources in Gaia. The dotted line
gives our colour limit, only stars to the left of which we consider. The solid
lines indicate the main-sequence region within which we consider HVS
candidates. The colour scale is logarithmic.

To provide perspective, we also overlay the density of HVS candidates in Gaia
in Fig. 14 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 anti-centre, 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
HEALPix222http://healpix.sf.net/ ($nside=3$; Górski et al., 2005; Zonca et
al., 2019) to divide the sky into equal area pixels. We determine the expected
number of HVSs using the simulations described in Section 4. We only observe
sources in HEALPix pixels where the ratio of the normalised number of expected
HVSs to candidates is larger than five333The HEALPix pixels we use are 0-19,
22-31, 35-42, 46-55, 58-62, 65, 66, 69-73, 78, 82-85, 93-96, 101, 102, 103 and
107 in the ring scheme based on RA and Dec..

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 centres444The LMC centre we use
has $l=280.4652$, $b=-32.8884$ deg. The SMC centre we define to be at
$l=302.8084$, $b=-44.3277$ deg..

The selection steps up to this point leave us with only 600 HVS candidates
with ${\rm G}<20$, 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.

Table 1: Overview of the HVS candidate selections used in this study. Selection | Section for reference
---|---
$\varpi_{\rm I}-\varpi_{\rm Gaia}<2\sigma$ | 2.2.1
$800<V_{\rm I}<3500$ $\text{km}\,\text{s}^{-1}$ | 2.2.1
$\varpi_{\rm I}/\sigma_{\varpi_{\rm I}}>5$ | 2.2.1
RUWE $<1.4$ | 2.2.1
$D_{\rm I}-D_{\rm BJ}<2\sigma$ | 2.2.1
$-1.5<M_{\rm G}-4.3\times({\rm G_{BP}-G_{RP}})<2.5$ | 2.2.2
$A_{\rm G}<1.5$ | 2.2.2
${\rm G_{BP}-G_{RP}}<0.5$ | 2.2.2
$N_{\rm HVS,sim}/N_{\rm HVS,candidates}>5$ | 2.2.3
Separation from LMC $>8$deg | 2.2.3
Separation from SMC $>3$deg | 2.2.3

The complete source list can be found
online555https://zenodo.org/doi/10.5281/zenodo.12179452. We provide the first
five sources in Table 2.

Table 2: Here we provide the first five sources, ordered by source_id, from our catalogue of 600 HVS candidates with their Gaia DR3 source_id, implied heliocentric radial velocity, and implied heliocentric distance. Gaia DR3 source_id |  $V_{\rm r,I}$ [$\text{km}\,\text{s}^{-1}$ ] | $D_{\rm I}[kpc]$
---|---|---
16647293239608960 | 805 | 6.6
57537894455885184 | 1429 | 15.2
292010016092347648 | 829 | 21.1
307220183209400448 | 638 | 12.2
331500110074917120 | 1092 | 18.7
… |  |

####  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-meter telescope located on the island of La Palma and the NTT is a
3.58-meter 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}<19$

  * •

$\textrm{Dec}>-20$ deg

  * •

$\textrm{RA}<60$ deg or $\textrm{RA}>230$ deg.

Our instrument specific selection for the NTT is

  * •

$\textrm{G}<19.3$

  * •

$\textrm{Dec}<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.

Table 3: Overview of the observing dates for our survey. Instrument | Observing dates | Run number
---|---|---
INT | 2022 Aug 20-21 | ING.NL.22B.002
INT | 2022 Sep 19-30 | ING.NL.22B.002
NTT | 2022 Nov 25-28 | 110.23SU.001 & 110.23SU.002
NTT | 2023 June 3 | 111.24MP.001
NTT | 2023 July 14-15 | 111.24MP.002

###  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
Å666https://www.ing.iac.es/astronomy/instruments/ids/idsgrat_tables.html. We
were allocated a total of twelve 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 setup gives us a resolution of about 500 at 3800
Å777https://www.eso.org/sci/facilities/lasilla/instruments/efosc/inst/Efosc2Grisms.html.
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<\lambda<5240$ Å, while for the INT we use the wavelength interval
$3600<\lambda<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%. We only
consider spectra with a 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 the appendix Table 5. We
will also consider these five sources in our observational sample, bringing
our total completeness to about 71%. In Fig. 3 we plot the sky distribution of
our HVS candidate catalogue, indicating which ones have been observed
successfully.

Figure 3: Sky distribution of our HVS candidate sample. Black points indicate
the stars with follow-up observations and the orange points the ones with
either SDSS or LAMOST radial velocity measurements. For the grey stars we do
not have radial velocity measurements.

We provide an observational log of our observations in an online
table888https://zenodo.org/doi/10.5281/zenodo.12179452, the first five rows of
which are shown in Table 4.

Table 4: Five sources from our observed HVS candidates in no particular order. The full table can be accessed online (see main text). The radial velocities are in the barycentric frame. Gaia DR3 source_id | radial velocity [$\text{km}\,\text{s}^{-1}$ ] | radial velocity uncertainty [$\text{km}\,\text{s}^{-1}$ ] | date of observation | instrument
---|---|---|---|---
6408940170344856320 | 220.6 | 38.0 | 2022-11-28 | NTT
16647293239608960 | 42.4 | 40.6 | 2022-09-23 | INT
4080886923180284928 | 32.8 | 35.4 | 2023-06-04 | NTT
3755045965083116160 | 281.0 | 21.2 | 2022-11-26 | NTT
4195180362312873984 | -75.7 | 80.2 | 2023-07-16 | NTT
… |  |  |  |

In Fig 4 we plot the implied radial velocity against the measured radial
velocity.

Figure 4: Implied radial velocity against the measured radial velocity for the
HVS candidates we observed. The dashed line gives the bisection. The blue
points indicate our measurements, the red star highlights our measurement of
S5-HVS1, and the orange points show the stars with radial velocity
measurements from SDSS and LAMOST.

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 can not 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 et al., 2009).
The velocity at which the HVS is ejected from the sphere of influence of the
MBH is given by

$V_{\text{ej}}=\sqrt{\frac{2Gm_{\text{c}}}{a}}\left(\frac{M_{\text{MBH}}}{M}\right)^{1/6},$
(8)

where we neglected a multiplicative constant factor of order of unity which
depends on the geometry of the encounter (see, e.g., fig.10 in Sari et al.
(2009). In eq. 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⊙, Eisenhauer et al., 2005; Ghez et al., 2008; Schödel et al., 2009;
Gillessen et al., 2009; Boehle et al., 2016; Gillessen et al., 2017; Do et
al., 2019; GRAVITY Collaboration et al., 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
Speedystar999https://github.com/fraserevans/speedystar (Contigiani et al.,
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 massive
black hole (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 characterised by the primary zero-age main sequence mass ($m_{\text{p}}$),
mass ratio ($q$), and the orbital semi-major 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⊙. The secondary is then
assigned a mass according to a mass ratio distribution of the form
$f(q)\propto q^{\gamma}$, with $0.1\leq q\leq 1$. Lastly, the orbital semi-
major 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 in Sgr A∗’s sphere of influence
(Bartko et al., 2010; Lu et al., 2013; Löckmann et al., 2010, 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 super-solar 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 et al. (2000) in AMUSE (Portegies Zwart
et al., 2009, 2013; Pelupessy 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 (Dotter, 2016; Choi et al., 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.

Firstly 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 et al., 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)101010https://github.com/gaia-unlimited/gaiaunlimited?tab=readme-ov-
file. Although Gaia is virtually complete down to our magnitude limit of ${\rm
G}=19.3$ except in crowded regions, we nonetheless 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 photo-
geometric 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 photo-geometric
distances for the simulated HVSs. After obtaining the photo-geometric
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 stars have
${\tt RUWE}>1.4$ (Penoyre et al., 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.

Figure 5: Observational completeness of our sample of HVS candidates as a
function of apparent G magnitude. We also show the distribution of simulated
HVSs and our observed candidates in G magnitude in red on the right axis.

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
prioritised 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.

Figure 6: Observational completeness of our sample of HVS candidates as a
function of the angle to the GC. In addition, we show the distribution of
simulated HVSs and our observed candidates as a function of their angular
separation from the GC in red on the right axis.

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}}<300$ $\text{km}\,\text{s}^{-1}$, which holds
for about 86% 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.

Figure 7: Predicted population of unbound HVSs in the Gaia catalogue. The
three columns display different assumptions on the MF of the primary of the
HVS progenitor binary population, as indicated above the respective columns.
The numbers per bin are calculated for a fiducial ejection rate of $10^{-6}$
yr-1. Top row: the distance-mass distribution for predicted HVSs in Gaia. The
black crosses indicate the median values of the respective panels, while the
grey crosses indicate those for the other panels for easy comparison. The
black star indicates the position of S5-HVS1. Bottom row: the extinction
corrected colour-magnitude distribution for predicted HVSs in Gaia.

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 on 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 over density of HVSs on the sky
is present. The cutoff 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 Appendix
Fig. 15, we additionally show the apparent magnitude distribution for a 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.

Figure 8: Completeness of unbound HVSs in the main Gaia catalogue as a
function of mass. The different lines correspond to different limits on the
heliocentric distance as indicated in the legend. Our observations are mostly
sensitive to HVSs within 50 kpc.

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 the completeness converges to nearly one, because those stars
do not live long enough to travel far from the GC. because of the effective
distance limit imposed by the short life times of massive stars.

###  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

$P(M|N_{\rm HVS})\propto\frac{\lambda^{N_{\rm HVS}}\exp{-\lambda}}{N_{\rm
HVS}!},$ (9)

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, marginalised
over the binary period and mass ratio distributions.

Figure 9: Posterior probability on the log of the HVS ejection rate and MF of
the primary in the progenitor binary. The solid line gives the $2\sigma$ upper
limit of the combined ejection rate and MF index. For comparison, the dotted
and dashed lines give the $2\sigma$ upper limit on the same parameters from
blind searches in the 6D phase-space Gaia eDR3 and DR3 samples respectively.
The eDR3 constraint was calculated in Evans et al. (2022b) and the DR3
constraint is calculated by multiplying the posteriors in Fig. 3 of Marchetti
et al. (2022) with the top-right panel of Fig. 5 in Evans et al. (2022b).

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 one 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. 16 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.

We can re-parameterise 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}]$

$N(M_{1}<M<M_{2})\propto\frac{1}{\kappa+1}\left(M_{2}^{\kappa+1}-M_{1}^{\kappa+1}\right).$
(10)

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.

Figure 10: Same as Fig. 9, but now only for HVSs with $M>1M_{\odot}$.

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.

Figure 11: Predicted distribution of flight times for HVSs to which our
observations are sensitive.

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.

###  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.

Figure 12: Number of unbound HVSs in the complete Gaia catalogue as a function
of the MF of the progenitor binary population.

Within $1\sigma$, we expect there to be between about 5 and 45 HVSs
($1.7<\kappa<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 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, secondly 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 parameterisations. 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 mass function. 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% 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 over-predict 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⊙,
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 et al., 2020; Shipp et
al., 2021; Koposov et al., 2023). In the run with the LMC, we found an
increase of about $2.8\%$ 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 S5 survey (Li et al., 2019a)
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 (Koposov
et al., 2020; Brown et al., 2014; Brown et al., 2018). An interesting question
to pose is then: what is the chance that the S5 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 S5 survey to be any stars with Gaia parallax
$\varpi<3\sigma_{\varpi}+0.2$, mock DECam (Dark Energy Survey Collaboration et
al., 2016) photometry $15<g<19.5$, $-0.4<(g-r)<0.1$, and falling within 2 deg
from an S5 pointing (see Li et al., 2019b, table 2) following Evans et al.
(2022b). Only HVSs with heliocentric radial velocities of $V_{\text{r}}>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 S5 is
about $0.14_{-0.11}^{+0.36}$. Although unlikely, the discovery of S5-HVS1 in
S5 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 to 5240
Å. We used this wavelength range in particular because the spectrum we
obtained for S5-HVS1 was much higher signal-to-noise 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⊙.
The difficulty with these stars is that due to their distance, their orbits
can not be confidently traced back to the GC (Irrgang et al., 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}>0$
deg for the Northern Galactic cap, and ${\rm Dec}>-15$ deg for the Southern
Galactic cap with a Galactic latitude $|l|>30$ deg (for the sky coverage see
Aihara et al., 2011). Brown et al. (2018) only studies 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% confidence
interval, we expect there to be between $0.3$ and $2.8$ HVSs in the MMT survey
footprint. For the 95% confidence interval, we expect there to be between
$0.1$ and $5.1$ HVSs respectively. Interestingly, there is thus a $>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 .

Figure 13: Flight time distribution for simulated HVSs matching the MMT survey
selections.

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\leq M_{\odot}\leq 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 a 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}}<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 recognised as such. If we additionally require an accurate
parallax measurement ($\varpi/\sigma_{\varpi}>5$), we only expect there to be
$1^{+1}_{-1}$ HVSs. If on top of that, we only consider HVSs with radial
velocities $>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}}<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.
Utilising 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% 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-100$ kpc.

## 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.

## Acknowledgements

The authors would like to thank C.A.L. Bailer Jones for kindly making the
source code to his work available.

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 et al., 2013, 2018,
2022), dustmaps (Green, 2018), healpy (Górski et al., 2005; Zonca et al.,
2019), sqlutilpy111111https://doi.org/10.5281/zenodo.6867957. This research
made use of ccdproc, an Astropy package for image reduction (Craig et al.,
2022).

## References

  * Abolfathi et al. (2018) Abolfathi B., et al., 2018, ApJS, 235, 42
  * Aihara et al. (2011) Aihara H., et al., 2011, ApJS, 193, 29
  * Asplund et al. (2009) Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, apj, 935, 167
  * Bailer-Jones (2015) Bailer-Jones C. A. L., 2015, PASP, 127, 994
  * Bailer-Jones et al. (2021) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147
  * Bartko et al. (2010) Bartko H., et al., 2010, ApJ, 708, 834
  * Boehle et al. (2016) Boehle A., et al., 2016, ApJ, 830, 17
  * Boubert et al. (2020) Boubert D., Erkal D., Gualandris A., 2020, MNRAS, 497, 2930
  * Brown (2015) Brown W. R., 2015, ARA&A, 53, 15
  * Brown et al. (2012) Brown W. R., Geller M. J., Kenyon S. J., 2012, ApJ, 751, 55
  * Brown et al. (2014) Brown W. R., Geller M. J., Kenyon S. J., 2014, ApJ, 787, 89
  * Brown et al. (2018) Brown W. R., Lattanzi M. G., Kenyon S. J., Geller M. J., 2018, ApJ, 866, 39
  * Cantat-Gaudin et al. (2023) Cantat-Gaudin T., et al., 2023, A&A, 669, A55
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Contigiani et al. (2019) Contigiani O., Rossi E. M., Marchetti T., 2019, MNRAS, 487, 4025
  * Cooper et al. (2023) Cooper A. P., et al., 2023, ApJ, 947, 37
  * Craig et al. (2022) Craig M., et al., 2022, astropy/ccdproc: 2.3.1 – fixes astropy 5.1 compatibility, doi:10.5281/zenodo.6533213, https://doi.org/10.5281/zenodo.6533213
  * Dalton et al. (2014) Dalton G., et al., 2014, in Ramsay S. K., McLean I. S., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V. p. 91470L (arXiv:1412.0843), doi:10.1117/12.2055132
  * Dark Energy Survey Collaboration et al. (2016) Dark Energy Survey Collaboration et al., 2016, MNRAS, 460, 1270
  * Do et al. (2015) Do T., Kerzendorf W., Winsor N., Støstad M., Morris M. R., Lu J. R., Ghez A. M., 2015, ApJ, 809, 143
  * Do et al. (2019) Do T., et al., 2019, Science, 365, 664
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Drimmel & Poggio (2018) Drimmel R., Poggio E., 2018, Research Notes of the AAS, 2, 210
  * Dunstall et al. (2015) Dunstall P. R., et al., 2015, A&A, 580, A93
  * Eisenhauer et al. (2005) Eisenhauer F., et al., 2005, The Astrophysical Journal, 628, 246
  * Eisenhauer et al. (2011) Eisenhauer F., et al., 2011, The Messenger, 143, 16
  * Erkal et al. (2019) Erkal D., et al., 2019, MNRAS, 487, 2685
  * Evans et al. (2022a) Evans F. A., Marchetti T., Rossi E. M., 2022a, MNRAS, 512, 2350
  * Evans et al. (2022b) Evans F. A., Marchetti T., Rossi E. M., 2022b, MNRAS, 517, 3469
  * Evans et al. (2023) Evans F. A., Rasskazov A., Remmelzwaal A., Marchetti T., Castro-Ginard A., Rossi E. M., Bovy J., 2023, MNRAS, 525, 561
  * Everall et al. (2021) Everall A., Boubert D., Koposov S. E., Smith L., Holl B., 2021, Monthly Notices of the Royal Astronomical Society, 502, 1908
  * Feldmeier-Krause et al. (2017) Feldmeier-Krause A., Kerzendorf W., Neumayer N., Schödel R., Nogueras-Lara F., Do T., de Zeeuw P. T., Kuntschner H., 2017, MNRAS, 464, 194
  * Fitzpatrick (1999) Fitzpatrick E. L., 1999, PASP, 111, 63
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration et al., 2018, A&A, 615, L15
  * GRAVITY Collaboration et al. (2019) GRAVITY Collaboration et al., 2019, A&A, 625, L10
  * Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
  * Gaia Collaboration et al. (2023) Gaia Collaboration et al., 2023, A&A, 674, A1
  * Generozov & Madigan (2020) Generozov A., Madigan A.-M., 2020, ApJ, 896, 137
  * Generozov & Perets (2022) Generozov A., Perets H. B., 2022, MNRAS, 513, 4257
  * Genzel et al. (2010) Genzel R., Eisenhauer F., Gillessen S., 2010, Reviews of Modern Physics, 82, 3121
  * Ghez et al. (2003) Ghez A. M., et al., 2003, The Astrophysical Journal, 586, L127
  * Ghez et al. (2008) Ghez A. M., et al., 2008, The Astrophysical Journal, 689, 1044
  * Gillessen et al. (2009) Gillessen S., Eisenhauer F., Trippe S., Alexander T., Genzel R., Martins F., Ott T., 2009, ApJ, 692, 1075
  * Gillessen et al. (2017) Gillessen S., et al., 2017, ApJ, 837, 30
  * Górski et al. (2005) Górski K. M., Hivon E., Banday A. J., Wandelt B. D., Hansen F. K., Reinecke M., Bartelmann M., 2005, ApJ, 622, 759
  * Gould & Quillen (2003) Gould A., Quillen A. C., 2003, The Astrophysical Journal, 592, 935
  * Green (2018) Green G., 2018, The Journal of Open Source Software, 3, 695
  * Guy et al. (2023) Guy J., et al., 2023, AJ, 165, 144
  * Habibi et al. (2017) Habibi M., et al., 2017, ApJ, 847, 120
  * Harris et al. (2020) Harris C. R., et al., 2020, Nature, 585, 357
  * Hills (1988) Hills J. G., 1988, Nature, 331, 687
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, Monthly Notices of the Royal Astronomical Society, 315, 543
  * Husser et al. (2013) Husser T. O., Wende-von Berg S., Dreizler S., Homeier D., Reiners A., Barman T., Hauschildt P. H., 2013, A&A, 553, A6
  * Irrgang et al. (2018) Irrgang A., Kreuzer S., Heber U., 2018, A&A, 620, A48
  * Irrgang et al. (2021) Irrgang A., Dimpel M., Heber U., Raddi R., 2021, A&A, 646, L4
  * Katz et al. (2023) Katz D., et al., 2023, A&A, 674, A5
  * Kenyon et al. (2018) Kenyon S. J., Bromley B. C., Brown W. R., Geller M. J., 2018, ApJ, 864, 130
  * Klessen et al. (2007) Klessen R. S., Spaans M., Jappsen A.-K., 2007, MNRAS, 374, L29
  * Kobayashi et al. (2012) Kobayashi S., Hainick Y., Sari R., Rossi E. M., 2012, The Astrophysical Journal, 748, 105
  * Koposov (2019) Koposov S. E., 2019, RVSpecFit: Radial velocity and stellar atmospheric parameter fitting, Astrophysics Source Code Library, record ascl:1907.013 (ascl:1907.013)
  * Koposov et al. (2020) Koposov S. E., et al., 2020, MNRAS, 491, 2465
  * Koposov et al. (2023) Koposov S. E., et al., 2023, MNRAS, 521, 4936
  * Kreuzer et al. (2020) Kreuzer S., Irrgang A., Heber U., 2020, A&A, 637, A53
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Li et al. (2019a) Li T. S., et al., 2019a, MNRAS, 490, 3508
  * Li et al. (2019b) Li T. S., et al., 2019b, MNRAS, 490, 3508
  * Löckmann et al. (2010) Löckmann U., Baumgardt H., Kroupa P., 2010, MNRAS, 402, 519
  * Lu et al. (2013) Lu J. R., Do T., Ghez A. M., Morris M. R., Yelda S., Matthews K., 2013, ApJ, 764, 155
  * Maness et al. (2007) Maness H., et al., 2007, ApJ, 669, 1024
  * Marchetti et al. (2018) Marchetti T., Contigiani O., Rossi E. M., Albert J. G., Brown A. G. A., Sesana A., 2018, MNRAS, 476, 4697
  * Marchetti et al. (2022) Marchetti T., Evans F. A., Rossi E. M., 2022, MNRAS, 515, 767
  * McMillan (2017) McMillan P. J., 2017, MNRAS, 465, 76
  * Moe & Stefano (2013) Moe M., Stefano R. D., 2013, The Astrophysical Journal, 778, 95
  * Moe & Stefano (2015) Moe M., Stefano R. D., 2015, The Astrophysical Journal, 810, 61
  * Nandakumar et al. (2018) Nandakumar G., Ryde N., Schultheis M., Thorsbro B., Jönsson H., Barklem P. S., Rich R. M., Fragkoudi F., 2018, MNRAS, 478, 4374
  * Nogueras-Lara et al. (2020) Nogueras-Lara F., et al., 2020, Nature Astronomy, 4, 377
  * Paumard et al. (2006) Paumard T., et al., 2006, ApJ, 643, 1011
  * Pelupessy et al. (2013) Pelupessy F. I., van Elteren, A. de Vries, N. McMillan, S. L. W. Drost, N. Portegies Zwart, S. F. 2013, A&A, 557, A84
  * Penoyre et al. (2022) Penoyre Z., Belokurov V., Evans N. W., 2022, MNRAS, 513, 2437
  * Portegies Zwart & McMillan (2018) Portegies Zwart S., McMillan S., 2018, Astrophysical Recipes; The art of AMUSE. IOP Publishing, doi:10.1088/978-0-7503-1320-9
  * Portegies Zwart et al. (2009) Portegies Zwart S., et al., 2009, New Astronomy, 14, 369
  * Portegies Zwart et al. (2013) Portegies Zwart S. F., McMillan S. L., van Elteren A., Pelupessy F. I., de Vries N., 2013, Computer Physics Communications, 184, 456
  * Przybilla et al. (2008) Przybilla N., Fernanda Nieva M., Heber U., Butler K., 2008, ApJ, 684, L103
  * Rossi et al. (2014) Rossi E. M., Kobayashi S., Sari R., 2014, The Astrophysical Journal, 795, 125
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sana et al. (2013) Sana H., et al., 2013, A&A, 550, A107
  * Sari et al. (2009) Sari R., Kobayashi S., Rossi E. M., 2009, The Astrophysical Journal, 708, 605
  * Schlafly & Finkbeiner (2011) Schlafly E. F., Finkbeiner D. P., 2011, The Astrophysical Journal, 737, 103
  * Schlegel et al. (1998) Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525
  * Schödel et al. (2009) Schödel R., Merritt D., Eckart A., 2009, A&A, 502, 91
  * Schultheis et al. (2019) Schultheis M., Rich R. M., Origlia L., Ryde N., Nandakumar G., Thorsbro B., Neumayer N., 2019, A&A, 627, A152
  * Shipp et al. (2021) Shipp N., et al., 2021, ApJ, 923, 149
  * Tody (1986) Tody D., 1986, in Crawford D. L., ed., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 627, Instrumentation in astronomy VI. p. 733, doi:10.1117/12.968154
  * Tody (1993) Tody D., 1993, in Hanisch R. J., Brissenden R. J. V., Barnes J., eds, Astronomical Society of the Pacific Conference Series Vol. 52, Astronomical Data Analysis Software and Systems II. p. 173
  * Vasiliev et al. (2020) Vasiliev E., Belokurov V., Erkal D., 2020, Monthly Notices of the Royal Astronomical Society, 501, 2279
  * Verberne et al. (2024) Verberne S., Koposov S. E., Rossi E. M., Marchetti T., Kuijken K., Penoyre Z., 2024, A&A, 684, A29
  * Virtanen et al. (2020) Virtanen P., et al., 2020, Nature Methods, 17, 261
  * Yu & Tremaine (2003) Yu Q., Tremaine S., 2003, ApJ, 599, 1129
  * Zhao et al. (2012) Zhao G., Zhao Y.-H., Chu Y.-Q., Jing Y.-P., Deng L.-C., 2012, Research in Astronomy and Astrophysics, 12, 723
  * Zonca et al. (2019) Zonca A., Singer L., Lenz D., Reinecke M., Rosset C., Hivon E., Gorski K., 2019, Journal of Open Source Software, 4, 1298
  * de Jong et al. (2019) de Jong R. S., et al., 2019, The Messenger, 175, 3
  * von Fellenberg et al. (2022) von Fellenberg S. D., et al., 2022, ApJ, 932, L6

##  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. 14 we show where HVS candidates in Gaia that match all but
our colour-magnitude selections end up in the HR diagram.

Figure 14: Same as Fig. 2, but now the overlain dashed line shows the density
of HVS candidates that meet all our selections, except the colour-magnitude
based ones.

##  Appendix B Radial velocities from SDSS and LAMOST

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

Table 5: Gaia DR3 source_id’s and radial velocities for the five stars for which we used literature radial velocity measurements. Gaia DR3 source_id | Radial velocity [$\text{km}\,\text{s}^{-1}$ ] | Source
---|---|---
2494761131458383872 | $-185\pm 13$ | SDSS
2700705779669486848 | $-73\pm 10$ | LAMOST
3792400463887610368 | $182\pm 13$ | LAMOST
4430178054799074688 | $197\pm 8$ | SDSS
4459380675615776640 | $-189\pm 5$ | SDSS

##  Appendix C G-magnitude distribution Gaia observable HVSs

In Fig. 15 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.

Figure 15: Apparent G-magnitude distribution of simulated HVSs predicted to be
observable by Gaia.

##  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. 16 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.

Figure 16: Posterior on the ejection rate and MF. Same as Fig. 9, but with a
log-uniform prior on the ejection rate of $1/\lambda$. The 95% upper limit is
indicated by the solid red line.