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: verberne@strw.leidenuniv.nl
(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 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. 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 MsubscriptMdirect-product\mathrm{M_{\odot}}roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT cannot be higher than 105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT yr-1 at 2σ2𝜎2\sigma2 italic_σ significance. Additionally, we predict that there are 5–45 unbound HVSs in the complete Gaia catalogue (1σ1𝜎1\sigma1 italic_σ 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: 2024pagerange: Hypervelocity star observations constrain the Galactic CentreD

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 10101010 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 R𝑅\vec{R}over→ start_ARG italic_R end_ARG be the vector from the GC to the star, and V𝑉\vec{V}over→ start_ARG italic_V end_ARG the star’s velocity vector relative to the GC. We can write these down in terms of the observables

  • R0subscript𝑅0\vec{R_{0}}over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG: Vector from the GC to the Sun

  • n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG: Unit vector from the Sun to the source

  • D𝐷Ditalic_D: Distance from the Sun to the source

  • V0subscript𝑉0\vec{V_{0}}over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG: Velocity of the Sun in the Galactic frame of rest

  • Vrsubscript𝑉rV_{\mathrm{r}}italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT: Radial velocity of the source relative to the Sun

  • μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG: Proper motion of the source in the Heliocentric frame (always perpendicular-to\perp to n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG)

giving

R=R0+Dn^𝑅subscript𝑅0𝐷^𝑛\vec{R}=\vec{R_{0}}+D\hat{n}over→ start_ARG italic_R end_ARG = over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG + italic_D over^ start_ARG italic_n end_ARG (1)

and

V=V0+Vrn^+Dμ.𝑉subscript𝑉0subscript𝑉r^𝑛𝐷𝜇\vec{V}=\vec{V_{0}}+V_{\rm r}\hat{n}+D\vec{\mu}.over→ start_ARG italic_V end_ARG = over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG + italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT over^ start_ARG italic_n end_ARG + italic_D over→ start_ARG italic_μ end_ARG . (2)

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

Refer to caption
Figure 1: A diagram of the position (R𝑅\vec{R}over→ start_ARG italic_R end_ARG) and velocity (V𝑉\vec{V}over→ start_ARG italic_V end_ARG), 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𝐷Ditalic_D) and its radial velocity (Vrsubscript𝑉𝑟V_{r}italic_V start_POSTSUBSCRIPT italic_r end_POSTSUBSCRIPT).

We use the measurement of R0subscript𝑅0\vec{R_{0}}over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG from GRAVITY Collaboration et al. (2018), which is 8.1228.1228.1228.122 kpc, and V0subscript𝑉0\vec{V_{0}}over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG from Drimmel & Poggio (2018), which is [12.9,245.6,7.78]12.9245.67.78[12.9,245.6,7.78][ 12.9 , 245.6 , 7.78 ] kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. The unit vector n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG is set by the sky coordinates of a star which are provided by Gaia, as is the proper motion vector μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG. For most sources in Gaia DR3, approximately 1.5 billion, we have measurements for all but D𝐷Ditalic_D and Vrsubscript𝑉rV_{\rm r}italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT. 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𝐷Ditalic_D and Vrsubscript𝑉rV_{\rm r}italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT 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 R𝑅\vec{R}over→ start_ARG italic_R end_ARG and velocity V𝑉\vec{V}over→ start_ARG italic_V end_ARG are parallel and thus obey

R×V=0.𝑅𝑉0\vec{R}\times\vec{V}=\vec{0}.over→ start_ARG italic_R end_ARG × over→ start_ARG italic_V end_ARG = over→ start_ARG 0 end_ARG . (3)

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

R0×V0+Vr(R0×n^)+D(R0×μ+n^×V0)+D2(n^×μ)=0subscript𝑅0subscript𝑉0subscript𝑉rsubscript𝑅0^𝑛𝐷subscript𝑅0𝜇^𝑛subscript𝑉0superscript𝐷2^𝑛𝜇0\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}over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG × over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG + italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT ( over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG × over^ start_ARG italic_n end_ARG ) + italic_D ( over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG × over→ start_ARG italic_μ end_ARG + over^ start_ARG italic_n end_ARG × over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) + italic_D start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( over^ start_ARG italic_n end_ARG × over→ start_ARG italic_μ end_ARG ) = over→ start_ARG 0 end_ARG (4)

and we can pull out the distance and radial velocity by taking the dot products with n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG and μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG and rearranging.

Taking the dot product of equation 4 with n^^𝑛\hat{n}over^ start_ARG italic_n end_ARG reveals

D=n^(R0×V0)R0(n^×μ)=V0μ(n^(R^0×V^0)R0^(n^×μ^))𝐷^𝑛subscript𝑅0subscript𝑉0subscript𝑅0^𝑛𝜇subscript𝑉0𝜇^𝑛subscript^𝑅0subscript^𝑉0^subscript𝑅0^𝑛^𝜇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)italic_D = divide start_ARG over^ start_ARG italic_n end_ARG ⋅ ( over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG × over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) end_ARG start_ARG over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ⋅ ( over^ start_ARG italic_n end_ARG × over→ start_ARG italic_μ end_ARG ) end_ARG = divide start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_μ end_ARG ( divide start_ARG over^ start_ARG italic_n end_ARG ⋅ ( over^ start_ARG italic_R end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT × over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG start_ARG over^ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ⋅ ( over^ start_ARG italic_n end_ARG × over^ start_ARG italic_μ end_ARG ) end_ARG ) (5)

and similarly the dot product of equation 4 with μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG gives

Vrsubscript𝑉r\displaystyle V_{{\rm r}}italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT =DV0(n^×μ)μ(R0×V0)R0(n^×μ)absent𝐷subscript𝑉0^𝑛𝜇𝜇subscript𝑅0subscript𝑉0subscript𝑅0^𝑛𝜇\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)}= divide start_ARG italic_D over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ⋅ ( over^ start_ARG italic_n end_ARG × over→ start_ARG italic_μ end_ARG ) - over→ start_ARG italic_μ end_ARG ⋅ ( over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG × over→ start_ARG italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ) end_ARG start_ARG over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ⋅ ( over^ start_ARG italic_n end_ARG × over→ start_ARG italic_μ end_ARG ) end_ARG (6)
=V0(DR0V^0(n^×μ^)μ^(R^0×V^0)R^0(n^×μ^)).absentsubscript𝑉0𝐷subscript𝑅0subscript^𝑉0^𝑛^𝜇^𝜇subscript^𝑅0subscript^𝑉0subscript^𝑅0^𝑛^𝜇\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).= italic_V start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( divide start_ARG divide start_ARG italic_D end_ARG start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⋅ ( over^ start_ARG italic_n end_ARG × over^ start_ARG italic_μ end_ARG ) - over^ start_ARG italic_μ end_ARG ⋅ ( over^ start_ARG italic_R end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT × over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG start_ARG over^ start_ARG italic_R end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⋅ ( over^ start_ARG italic_n end_ARG × over^ start_ARG italic_μ end_ARG ) end_ARG ) .

This means that for a star moving on a radial trajectory, the distance and radial velocity are set by the observables available for 1.5similar-toabsent1.5\sim 1.5∼ 1.5B 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, DIsubscript𝐷ID_{\rm I}italic_D start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT, and implied radial velocity, Vr,Isubscript𝑉rIV_{{\rm r,I}}italic_V start_POSTSUBSCRIPT roman_r , roman_I end_POSTSUBSCRIPT.

We convert DIsubscript𝐷ID_{\rm I}italic_D start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT to implied parallax (ϖIsubscriptitalic-ϖI\varpi_{\rm I}italic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT) to allow for convenient comparison to measurements by taking the inverse of DIsubscript𝐷ID_{\rm I}italic_D start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT. We calculate the uncertainty on ϖIsubscriptitalic-ϖI\varpi_{\rm I}italic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT 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 ϖμproportional-toitalic-ϖ𝜇\varpi\propto\muitalic_ϖ ∝ italic_μ, we can approximate the uncertainty on ϖitalic-ϖ\varpiitalic_ϖ as

σϖ=[R0(n^×μ^δ)σμδ]2+[R0(n^×μ^α)σμα]2n^(R0×V^0),subscript𝜎italic-ϖsuperscriptdelimited-[]subscript𝑅0^𝑛subscript^𝜇𝛿subscript𝜎subscript𝜇𝛿2superscriptdelimited-[]subscript𝑅0^𝑛subscript^𝜇𝛼subscript𝜎subscript𝜇𝛼2^𝑛subscript𝑅0subscript^𝑉0\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)},italic_σ start_POSTSUBSCRIPT italic_ϖ end_POSTSUBSCRIPT = divide start_ARG square-root start_ARG [ over→ start_ARG italic_R end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ⋅ ( over^ start_ARG italic_n end_ARG × over^ start_ARG italic_μ end_ARG start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT ) italic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT + [ over→ start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG ⋅ ( over^ start_ARG italic_n end_ARG × over^ start_ARG italic_μ end_ARG start_POSTSUBSCRIPT italic_α ∗ end_POSTSUBSCRIPT ) italic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_α ∗ end_POSTSUBSCRIPT end_POSTSUBSCRIPT ] start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG start_ARG over^ start_ARG italic_n end_ARG ⋅ ( over→ start_ARG italic_R end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT × over^ start_ARG italic_V end_ARG start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ) end_ARG , (7)

with μ^δ\hat{\mu}{{}_{\delta}}over^ start_ARG italic_μ end_ARG start_FLOATSUBSCRIPT italic_δ end_FLOATSUBSCRIPT and μ^αsubscript^𝜇𝛼\hat{\mu}_{\alpha*}over^ start_ARG italic_μ end_ARG start_POSTSUBSCRIPT italic_α ∗ end_POSTSUBSCRIPT the basis vectors for the proper motion in declination (Dec) and right ascension (RA) respectively, and σμδsubscript𝜎subscript𝜇𝛿\sigma_{\mu_{\delta}}italic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_δ end_POSTSUBSCRIPT end_POSTSUBSCRIPT and σμαsubscript𝜎subscript𝜇𝛼\sigma_{\mu_{\alpha*}}italic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_α ∗ end_POSTSUBSCRIPT end_POSTSUBSCRIPT the uncertainties on the proper motions in Dec and RA respectively (where we take σμα=σμαcos2δsubscript𝜎subscript𝜇𝛼subscript𝜎subscript𝜇𝛼superscript2𝛿\sigma_{\mu_{\alpha}*}=\sigma_{\mu_{\alpha}}\cos^{2}\deltaitalic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT ∗ end_POSTSUBSCRIPT = italic_σ start_POSTSUBSCRIPT italic_μ start_POSTSUBSCRIPT italic_α end_POSTSUBSCRIPT end_POSTSUBSCRIPT roman_cos start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_δ). For a more accurate uncertainty analysis one would sample over the uncertainty on all observables.

2.2 HVS candidate selection criteria

The ϖIsubscriptitalic-ϖI\varpi_{\rm I}italic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT and Vr,Isubscript𝑉rIV_{\rm r,I}italic_V start_POSTSUBSCRIPT roman_r , roman_I end_POSTSUBSCRIPT 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, VIsubscript𝑉IV_{\text{I}}italic_V start_POSTSUBSCRIPT I end_POSTSUBSCRIPT, in the range [800,3500]8003500[800,3500][ 800 , 3500 ] kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT. 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σ2𝜎2\sigma2 italic_σ, 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 ϖIsubscriptitalic-ϖI\varpi_{\text{I}}italic_ϖ start_POSTSUBSCRIPT I end_POSTSUBSCRIPT is the measured uncertainty on μ𝜇\vec{\mu}over→ start_ARG italic_μ end_ARG from Gaia. Because comparing parallaxes is not constraining when both are highly uncertain, we only consider sources where ϖIsubscriptitalic-ϖI\varpi_{\text{I}}italic_ϖ start_POSTSUBSCRIPT I end_POSTSUBSCRIPT 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σ2𝜎2\sigma2 italic_σ, where we use the 16th and 84th percentiles from Bailer-Jones et al. (2021) as the negative and positive 1σ1𝜎1\sigma1 italic_σ 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, GRPsubscriptGRP{\rm G_{RP}}roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT, and GBPsubscriptGBP{\rm G_{BP}}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT 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 RV=AV/E(BV)subscript𝑅Vsubscript𝐴VEBVR_{\rm V}=A_{\rm V}/{\rm E(B-V)}italic_R start_POSTSUBSCRIPT roman_V end_POSTSUBSCRIPT = italic_A start_POSTSUBSCRIPT roman_V end_POSTSUBSCRIPT / roman_E ( roman_B - roman_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 ΔMG=2ΔsubscriptMG2\Delta{\rm M_{G}}=2roman_Δ roman_M start_POSTSUBSCRIPT roman_G end_POSTSUBSCRIPT = 2 of the main-sequence, which we define to be MG=4.3×(GBPGRP)+0.5subscriptMG4.3subscriptGBPsubscriptGRP0.5{\rm M_{G}}=4.3\times({\rm G_{BP}-G_{RP}})+0.5roman_M start_POSTSUBSCRIPT roman_G end_POSTSUBSCRIPT = 4.3 × ( roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT ) + 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 AG<1.5subscript𝐴G1.5A_{\text{G}}<1.5italic_A start_POSTSUBSCRIPT G end_POSTSUBSCRIPT < 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 GBPGRPsubscriptGBPsubscriptGRP{\rm G_{BP}-G_{RP}}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT colour below 0.5, which corresponds to a mass 1.3greater-than-or-equivalent-toabsent1.3\gtrsim 1.3≳ 1.3 M. We show these selections in Fig. 2 in combination with a reference HR diagram for sources in Gaia.

Refer to caption
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𝑛𝑠𝑖𝑑𝑒3nside=3italic_n italic_s italic_i italic_d italic_e = 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𝑙280.4652l=280.4652italic_l = 280.4652, b=32.8884𝑏32.8884b=-32.8884italic_b = - 32.8884 deg. The SMC centre we define to be at l=302.8084𝑙302.8084l=302.8084italic_l = 302.8084, b=44.3277𝑏44.3277b=-44.3277italic_b = - 44.3277 deg..

The selection steps up to this point leave us with only 600 HVS candidates with G<20G20{\rm G}<20roman_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
ϖIϖGaia<2σsubscriptitalic-ϖIsubscriptitalic-ϖGaia2𝜎\varpi_{\rm I}-\varpi_{\rm Gaia}<2\sigmaitalic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT - italic_ϖ start_POSTSUBSCRIPT roman_Gaia end_POSTSUBSCRIPT < 2 italic_σ 2.2.1
800<VI<3500800subscript𝑉I3500800<V_{\rm I}<3500800 < italic_V start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT < 3500 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 2.2.1
ϖI/σϖI>5subscriptitalic-ϖIsubscript𝜎subscriptitalic-ϖI5\varpi_{\rm I}/\sigma_{\varpi_{\rm I}}>5italic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT / italic_σ start_POSTSUBSCRIPT italic_ϖ start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT end_POSTSUBSCRIPT > 5 2.2.1
RUWE <1.4absent1.4<1.4< 1.4 2.2.1
DIDBJ<2σsubscript𝐷Isubscript𝐷BJ2𝜎D_{\rm I}-D_{\rm BJ}<2\sigmaitalic_D start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT - italic_D start_POSTSUBSCRIPT roman_BJ end_POSTSUBSCRIPT < 2 italic_σ 2.2.1
1.5<MG4.3×(GBPGRP)<2.51.5subscript𝑀G4.3subscriptGBPsubscriptGRP2.5-1.5<M_{\rm G}-4.3\times({\rm G_{BP}-G_{RP}})<2.5- 1.5 < italic_M start_POSTSUBSCRIPT roman_G end_POSTSUBSCRIPT - 4.3 × ( roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT ) < 2.5 2.2.2
AG<1.5subscript𝐴G1.5A_{\rm G}<1.5italic_A start_POSTSUBSCRIPT roman_G end_POSTSUBSCRIPT < 1.5 2.2.2
GBPGRP<0.5subscriptGBPsubscriptGRP0.5{\rm G_{BP}-G_{RP}}<0.5roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT < 0.5 2.2.2
NHVS,sim/NHVS,candidates>5subscript𝑁HVSsimsubscript𝑁HVScandidates5N_{\rm HVS,sim}/N_{\rm HVS,candidates}>5italic_N start_POSTSUBSCRIPT roman_HVS , roman_sim end_POSTSUBSCRIPT / italic_N start_POSTSUBSCRIPT roman_HVS , roman_candidates end_POSTSUBSCRIPT > 5 2.2.3
Separation from LMC >8absent8>8> 8deg 2.2.3
Separation from SMC >3absent3>3> 3deg 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 Vr,Isubscript𝑉rIV_{\rm r,I}italic_V start_POSTSUBSCRIPT roman_r , roman_I end_POSTSUBSCRIPT [kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ] DI[kpc]subscript𝐷Idelimited-[]𝑘𝑝𝑐D_{\rm I}[kpc]italic_D start_POSTSUBSCRIPT roman_I end_POSTSUBSCRIPT [ italic_k italic_p italic_c ]
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

  • G<19G19\textrm{G}<19G < 19

  • Dec>20Dec20\textrm{Dec}>-20Dec > - 20 deg

  • RA<60RA60\textrm{RA}<60RA < 60 deg or RA>230RA230\textrm{RA}>230RA > 230 deg.

Our instrument specific selection for the NTT is

  • G<19.3G19.3\textrm{G}<19.3G < 19.3

  • Dec<0Dec0\textrm{Dec}<0Dec < 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 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT.

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 R=λΔλ=1600R𝜆Δ𝜆1600\textrm{R}=\frac{\lambda}{\Delta\lambda}=1600R = divide start_ARG italic_λ end_ARG start_ARG roman_Δ italic_λ end_ARG = 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<λ<52403500𝜆52403500<\lambda<52403500 < italic_λ < 5240 Å, while for the INT we use the wavelength interval 3600<λ<68003600𝜆68003600<\lambda<68003600 < italic_λ < 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.

Refer to caption
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 [kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ] radial velocity uncertainty [kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ] 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.

Refer to caption
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

Vej=2Gmca(MMBHM)1/6,subscript𝑉ej2𝐺subscript𝑚c𝑎superscriptsubscript𝑀MBH𝑀16V_{\text{ej}}=\sqrt{\frac{2Gm_{\text{c}}}{a}}\left(\frac{M_{\text{MBH}}}{M}% \right)^{1/6},italic_V start_POSTSUBSCRIPT ej end_POSTSUBSCRIPT = square-root start_ARG divide start_ARG 2 italic_G italic_m start_POSTSUBSCRIPT c end_POSTSUBSCRIPT end_ARG start_ARG italic_a end_ARG end_ARG ( divide start_ARG italic_M start_POSTSUBSCRIPT MBH end_POSTSUBSCRIPT end_ARG start_ARG italic_M end_ARG ) start_POSTSUPERSCRIPT 1 / 6 end_POSTSUPERSCRIPT , (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, mcsubscript𝑚cm_{\text{c}}italic_m start_POSTSUBSCRIPT c end_POSTSUBSCRIPT is the mass of the captured companion, MMBHsubscript𝑀MBHM_{\text{MBH}}italic_M start_POSTSUBSCRIPT MBH end_POSTSUBSCRIPT in the case of the Galaxy is the mass of Sgr A (4×1064superscript1064\times 10^{6}4 × 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT 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𝑀Mitalic_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 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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 (mpsubscript𝑚pm_{\text{p}}italic_m start_POSTSUBSCRIPT p end_POSTSUBSCRIPT), mass ratio (q𝑞qitalic_q), and the orbital semi-major axis (a𝑎aitalic_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(mp)mpκproportional-to𝑓subscript𝑚psuperscriptsubscript𝑚p𝜅f(m_{\text{p}})\propto m_{\text{p}}^{-\kappa}italic_f ( italic_m start_POSTSUBSCRIPT p end_POSTSUBSCRIPT ) ∝ italic_m start_POSTSUBSCRIPT p end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - italic_κ end_POSTSUPERSCRIPT in the range [0.1,100]0.1100[0.1,100][ 0.1 , 100 ] M. The secondary is then assigned a mass according to a mass ratio distribution of the form f(q)qγproportional-to𝑓𝑞superscript𝑞𝛾f(q)\propto q^{\gamma}italic_f ( italic_q ) ∝ italic_q start_POSTSUPERSCRIPT italic_γ end_POSTSUPERSCRIPT, with 0.1q10.1𝑞10.1\leq q\leq 10.1 ≤ italic_q ≤ 1. Lastly, the orbital semi-major axis are assigned assuming a binary orbital period distribution of the form f(logP)(logP)πproportional-to𝑓𝑃superscript𝑃𝜋f(\log P)\propto(\log P)^{\pi}italic_f ( roman_log italic_P ) ∝ ( roman_log italic_P ) start_POSTSUPERSCRIPT italic_π end_POSTSUPERSCRIPT. Different values for the three parameters κ𝜅\kappaitalic_κ, γ𝛾\gammaitalic_γ, and π𝜋\piitalic_π are explored with uniform priors to account for uncertainty in the HVS progenitor population. For κ𝜅\kappaitalic_κ we use the interval [0.3,3.3]0.33.3[0.3,3.3][ 0.3 , 3.3 ] to allow for both very top-heavy and bottom-heavy MFs. For γ𝛾\gammaitalic_γ we use the interval [2,2]22[-2,2][ - 2 , 2 ] and for π𝜋\piitalic_π the interval [0,2]02[0,2][ 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 π𝜋\piitalic_π 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𝑞qitalic_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=0.0142subscript𝑍direct-product0.0142Z_{\odot}=0.0142italic_Z start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT = 0.0142 (Asplund et al., 2009). Additionally, we performed tests with super-solar metallicity values of [Fe/H]=0.3delimited-[]FeH0.3[\mathrm{Fe/H}]=0.3[ roman_Fe / roman_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 G=20.7G20.7\text{G}=20.7G = 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 G=19.3G19.3{\rm G}=19.3roman_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 𝚁𝚄𝚆𝙴>1.4𝚁𝚄𝚆𝙴1.4{\tt RUWE}>1.4typewriter_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.

Refer to caption
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.

Refer to caption
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 G=0G0{\rm G}=0roman_G = 0 and G=19.4G19.4{\rm G}=19.4roman_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 Vr,IVr<300subscript𝑉rIsubscript𝑉r300V_{{\rm r,I}}-V_{{\rm r}}<300italic_V start_POSTSUBSCRIPT roman_r , roman_I end_POSTSUBSCRIPT - italic_V start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT < 300 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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.720.720.720.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.

Refer to caption
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 106superscript10610^{-6}10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 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, κ=1.7𝜅1.7\kappa=1.7italic_κ = 1.7, and one close to the canonical Salpeter MF, κ=2.3𝜅2.3\kappa=2.3italic_κ = 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 5similar-toabsent5\sim 5∼ 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 8similar-toabsent8\sim 8∼ 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 κ=2.3𝜅2.3\kappa=2.3italic_κ = 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.

Refer to caption
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 1Msimilar-toabsent1subscriptMdirect-product\sim 1\;{\rm M_{\odot}}∼ 1 roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT 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𝑀Mitalic_M is therefore given by

P(M|NHVS)λNHVSexpλNHVS!,proportional-to𝑃conditional𝑀subscript𝑁HVSsuperscript𝜆subscript𝑁HVS𝜆subscript𝑁HVSP(M|N_{\rm HVS})\propto\frac{\lambda^{N_{\rm HVS}}\exp{-\lambda}}{N_{\rm HVS}!},italic_P ( italic_M | italic_N start_POSTSUBSCRIPT roman_HVS end_POSTSUBSCRIPT ) ∝ divide start_ARG italic_λ start_POSTSUPERSCRIPT italic_N start_POSTSUBSCRIPT roman_HVS end_POSTSUBSCRIPT end_POSTSUPERSCRIPT roman_exp - italic_λ end_ARG start_ARG italic_N start_POSTSUBSCRIPT roman_HVS end_POSTSUBSCRIPT ! end_ARG , (9)

with NHVSsubscript𝑁HVSN_{\rm HVS}italic_N start_POSTSUBSCRIPT roman_HVS end_POSTSUBSCRIPT the number of observed HVSs and λ𝜆\lambdaitalic_λ the average number of expected HVSs according to the simulations for a given set of parameters M𝑀Mitalic_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.

Refer to caption
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σ2𝜎2\sigma2 italic_σ upper limit of the combined ejection rate and MF index. For comparison, the dotted and dashed lines give the 2σ2𝜎2\sigma2 italic_σ 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 34similar-toabsent34\sim 34∼ 34M 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/λ1𝜆1/\lambda1 / italic_λ. The posterior with the log-uniform prior shows that the upper limit remains, while allowing lower ejection rates with equal probability.

For κ1.5greater-than-or-equivalent-to𝜅1.5\kappa\gtrsim 1.5italic_κ ≳ 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]M0.1100subscriptMdirect-product[0.1,100]\ \mathrm{M}_{\odot}[ 0.1 , 100 ] roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT, 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 [M1,M2]subscript𝑀1subscript𝑀2[M_{1},M_{2}][ italic_M start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_M start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ]

N(M1<M<M2)1κ+1(M2κ+1M1κ+1).proportional-to𝑁subscript𝑀1𝑀subscript𝑀21𝜅1superscriptsubscript𝑀2𝜅1superscriptsubscript𝑀1𝜅1N(M_{1}<M<M_{2})\propto\frac{1}{\kappa+1}\left(M_{2}^{\kappa+1}-M_{1}^{\kappa+% 1}\right).italic_N ( italic_M start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT < italic_M < italic_M start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ) ∝ divide start_ARG 1 end_ARG start_ARG italic_κ + 1 end_ARG ( italic_M start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_κ + 1 end_POSTSUPERSCRIPT - italic_M start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_κ + 1 end_POSTSUPERSCRIPT ) . (10)

The ejection rate for stars in the mass range [M1,M2]subscript𝑀1subscript𝑀2[M_{1},M_{2}][ italic_M start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_M start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ] is then the total ejection rate times the fraction of stars in the mass interval [M1,M2]subscript𝑀1subscript𝑀2[M_{1},M_{2}][ italic_M start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT , italic_M start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT ]. In Fig. 10 we give the resulting posterior for HVSs more massive than 1M1subscriptMdirect-product1\mathrm{M}_{\odot}1 roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT, making them potentially accessible to Gaia.

Refer to caption
Figure 10: Same as Fig. 9, but now only for HVSs with M>1M𝑀1subscript𝑀direct-productM>1M_{\odot}italic_M > 1 italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT.

The figure demonstrates that irrespective of the MF, the ejection rate of HVSs more massive than 1M1subscriptMdirect-product1\mathrm{M}_{\odot}1 roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT cannot be higher than about 105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT yr-1 at 2σ2𝜎2\sigma2 italic_σ confidence.

In Fig. 11 we show the flight time (i.e. time since disruption) distribution of HVSs that our observations are sensitive to.

Refer to caption
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 50100similar-toabsent50100\sim 50-100∼ 50 - 100 Myr. Our ejection rate constraint therefore applies to ejections within the last 50100similar-toabsent50100\sim 50-100∼ 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.

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

Within 1σ1𝜎1\sigma1 italic_σ, we expect there to be between about 5 and 45 HVSs (1.7<κ<2.31.7𝜅2.31.7<\kappa<2.31.7 < italic_κ < 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 α𝛼\alphaitalic_α (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 [Fe/H]=0.3delimited-[]FeH0.3[{\rm Fe/H}]=0.3[ roman_Fe / roman_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.010.10.010.10.01-0.10.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 κ=2.3𝜅2.3\kappa=2.3italic_κ = 2.3 including the LMC with a mass of 1.5×10111.5superscript10111.5\times 10^{11}1.5 × 10 start_POSTSUPERSCRIPT 11 end_POSTSUPERSCRIPT 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%percent2.82.8\%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 ϖ<3σϖ+0.2italic-ϖ3subscript𝜎italic-ϖ0.2\varpi<3\sigma_{\varpi}+0.2italic_ϖ < 3 italic_σ start_POSTSUBSCRIPT italic_ϖ end_POSTSUBSCRIPT + 0.2, mock DECam (Dark Energy Survey Collaboration et al., 2016) photometry 15<g<19.515𝑔19.515<g<19.515 < italic_g < 19.5, 0.4<(gr)<0.10.4𝑔𝑟0.1-0.4<(g-r)<0.1- 0.4 < ( italic_g - italic_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 Vr>800subscript𝑉r800V_{\text{r}}>800italic_V start_POSTSUBSCRIPT r end_POSTSUBSCRIPT > 800 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT are considered, since these were inspected in Koposov et al. (2020). We find that within the 1σ1𝜎1\sigma1 italic_σ posterior probability on the MF and ejection rate from Fig. 9, the number of expected HVSs in S5 is about 0.140.11+0.36superscriptsubscript0.140.110.360.14_{-0.11}^{+0.36}0.14 start_POSTSUBSCRIPT - 0.11 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT + 0.36 end_POSTSUPERSCRIPT. 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±12plus-or-minus99512995\pm 12995 ± 12 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT 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±2.7plus-or-minus10172.71017\pm 2.71017 ± 2.7 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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 Dec>0Dec0{\rm Dec}>0roman_Dec > 0 deg for the Northern Galactic cap, and Dec>15Dec15{\rm Dec}>-15roman_Dec > - 15 deg for the Southern Galactic cap with a Galactic latitude |l|>30𝑙30|l|>30| italic_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 275275275275 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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]1.72.3[1.7,2.3][ 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.30.30.30.3 and 2.82.82.82.8 HVSs in the MMT survey footprint. For the 95% confidence interval, we expect there to be between 0.10.10.10.1 and 5.15.15.15.1 HVSs respectively. Interestingly, there is thus a >2σabsent2𝜎>2\sigma> 2 italic_σ 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 .

Refer to caption
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.5M42.5subscript𝑀direct-product42.5\leq M_{\odot}\leq 42.5 ≤ italic_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT ≤ 4 to be 1.5×1061.5superscript1061.5\times 10^{-6}1.5 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT 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σ2𝜎2\sigma2 italic_σ upper limit peaks at about 1.3×1061.3superscript1061.3\times 10^{-6}1.3 × 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT yr-1 for a MF index κ2.5similar-to𝜅2.5\kappa\sim 2.5italic_κ ∼ 2.5, while the mode of the posterior peaks at about 2.8×1072.8superscript1072.8\times 10^{-7}2.8 × 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT 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<RVS16.2{}_{\mathrm{RVS}}<16.2start_FLOATSUBSCRIPT roman_RVS end_FLOATSUBSCRIPT < 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 32+4subscriptsuperscript3423^{+4}_{-2}3 start_POSTSUPERSCRIPT + 4 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 2 end_POSTSUBSCRIPT HVSs within the Gaia DR4 radial velocity catalogue if the MF index is between 1.71.71.71.7 and 2.32.32.32.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 (ϖ/σϖ>5italic-ϖsubscript𝜎italic-ϖ5\varpi/\sigma_{\varpi}>5italic_ϖ / italic_σ start_POSTSUBSCRIPT italic_ϖ end_POSTSUBSCRIPT > 5), we only expect there to be 11+1subscriptsuperscript1111^{+1}_{-1}1 start_POSTSUPERSCRIPT + 1 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 1 end_POSTSUBSCRIPT HVSs. If on top of that, we only consider HVSs with radial velocities >500absent500>500> 500 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, 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 GBP/GRPsubscriptGBPsubscriptGRP{\rm G_{BP}/G_{RP}}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT / roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT 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 (GBPGRP0.7less-than-or-similar-tosubscriptGBPsubscriptGRP0.7{\rm G_{BP}-G_{RP}\lesssim 0.7}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT ≲ 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 GBP/GRPsubscriptGBPsubscriptGRP{\rm G_{BP}/G_{RP}}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT / roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT spectra, we only consider HVSs with radial velocities greater than 1000 kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT, in addition to only selecting stars with GBPGRP<0.7subscriptGBPsubscriptGRP0.7{\rm G_{BP}-G_{RP}}<0.7roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT - roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT < 0.7. Within this parameter range, we expect there to be 75+9subscriptsuperscript7957^{+9}_{-5}7 start_POSTSUPERSCRIPT + 9 end_POSTSUPERSCRIPT start_POSTSUBSCRIPT - 5 end_POSTSUBSCRIPT HVSs in Gaia DR4. However, significant improvements would be required in the technique of obtaining these radial velocities from the GBP/GRPsubscriptGBPsubscriptGRP{\rm G_{BP}/G_{RP}}roman_G start_POSTSUBSCRIPT roman_BP end_POSTSUBSCRIPT / roman_G start_POSTSUBSCRIPT roman_RP end_POSTSUBSCRIPT 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 0002000020\,00020 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 1M1subscriptMdirect-product1\mathrm{M_{\odot}}1 roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT. We are therefore able to provide robust constraints on the ejection rate for HVSs more massive than 1M1subscriptMdirect-product1\mathrm{M_{\odot}}1 roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT independent on the MF. Over the past 50100similar-toabsent50100\sim 50-100∼ 50 - 100 Myr, the average ejection rate for these stars cannot have been higher than about 105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT yr-1 at 2σ2𝜎2\sigma2 italic_σ 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 150similar-toabsent150\sim 150∼ 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 MsubscriptMdirect-product{\rm M_{\odot}}roman_M start_POSTSUBSCRIPT ⊙ end_POSTSUBSCRIPT and be at Heliocentric distances between 101001010010-10010 - 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.

Refer to caption
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 [kms1kmsuperscripts1\text{km}\,\text{s}^{-1}km s start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ] Source
2494761131458383872 185±13plus-or-minus18513-185\pm 13- 185 ± 13 SDSS
2700705779669486848 73±10plus-or-minus7310-73\pm 10- 73 ± 10 LAMOST
3792400463887610368 182±13plus-or-minus18213182\pm 13182 ± 13 LAMOST
4430178054799074688 197±8plus-or-minus1978197\pm 8197 ± 8 SDSS
4459380675615776640 189±5plus-or-minus1895-189\pm 5- 189 ± 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 κ=2.3𝜅2.3\kappa=2.3italic_κ = 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.

Refer to caption
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 λ𝜆\lambdaitalic_λ. In Fig. 16 we instead show the posterior if we use a log-uniform prior equal to 1/λ1𝜆1/\lambda1 / italic_λ. 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.

Refer to caption
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/λ1𝜆1/\lambda1 / italic_λ. The 95% upper limit is indicated by the solid red line.