© The Authors 2026

1 Introduction

In the Galactic stellar halo, most stars exhibit random motions with velocities of 100-150 km s⁻¹. However, some high-velocity stars travel significantly faster - on the order of several hundred kilometers per second - with a subset known as hypervelocity stars (HVSs) even exceeding the local escape speed. Such objects are indicative of extreme dynamical processes. The classic Hills mechanism (Hills 1988) was the first to predict the existence of HVSs. In this scenario, when a stellar binary passes sufficiently close to the massive black hole (MBH) at the Galactic center (GC), tidal forces disrupt the binary. Ultimately, one component is captured by the MBH, while its companion is ejected at speeds potentially exceeding 1000 km s⁻¹.

Subsequent work has shown that single stars can also achieve extremely high velocities through interactions with more complex MBH systems. For example, encounters between a star and a binary MBH or an MBH-intermediate-mass black hole (IMBH) pair can efficiently transfer energy to the star, ejecting it at hundreds to thousands of kilometers per second (Yu & Tremaine 2003; Gualandris & Merritt 2009; Rasskazov et al. 2019). Similarly, a dense cluster of stellar-mass black holes orbiting the central MBH can scatter passing stars so that they become HVSs (O’Leary & Loeb 2008). Infalling globular clusters, when tidally disrupted by the GC’s MBH (or an MBH binary), can also eject stars at extreme velocities (Capuzzo-Dolcetta & Fragione 2015; Fragione et al. 2017).

Supernovae in binary systems constitute another natural channel for producing fast runaway stars. Core-collapse supernovae disrupt their progenitor binaries and can impart kicks of up to 300-400 km s⁻¹ (Evans et al. 2020), accounting for many OB runaway stars found above the Galactic plane (Blaauw 1961; Portegies Zwart 2000). Thermonuclear supernovae, in contrast, can accelerate companions to much higher speeds. In particular, detonations in white dwarf-white dwarf binaries or in systems pairing a white dwarf with a nondegenerate star have been proposed to yield velocities exceeding 1000 km s⁻¹ (Shen et al. 2018; Wang & Han 2009; Bauer et al. 2019; Neunteufel 2020; El-Badry et al. 2023).

Dense stellar environments and galactic interactions provide alternative channels for hypervelocity ejections. Close multibody encounters within young massive clusters can scatter stars onto unbound trajectories (Perets & Subr 2012; Oh & Kroupa 2016). Specifically, three-body interactions - particularly single star-binary encounters involving compact objects - can eject stars at velocities reaching several hundred kilometers per second (Cabrera & Rodriguez 2023; Weatherford et al. 2023; Grondin et al. 2023, 2024; Evans et al. 2025). Additionally, the tidal disruption of infalling dwarf galaxies can strip and accelerate stars to extreme velocities (Abadi et al. 2009; Piffl et al. 2011). Collectively, these mechanisms illustrate the diverse dynamical processes capable of populating the Galactic halo with HVSs.

The first candidate HVS was identified by Brown et al. (2005). It was a B-type star traveling at a total velocity of 709 km s⁻¹ through the Galactic halo. After that, over two dozen HVSs, all early-type, have been identified either by incidental discoveries (Hirsch et al. 2005; Heber et al. 2008; Koposov et al. 2020) or follow-up surveys (Brown et al. 2014). More recently, the Southern Stellar Stream Spectroscopic Survey (S5) reported the fastest known HVS to date (Koposov et al. 2020). This A-type star moves at 1755 km s⁻¹ and lies at a heliocentric distance of approximately 9 kpc, with a GC-related origin according to orbital tracing.

The advent of the Gaia mission has revolutionized the search for HVSs by providing accurate astrometry. Thanks to Gaia Data Releases 2 and 3 (Gaia Collaboration 2018a, 2021), dozens of new HVS candidates have been uncovered (e.g., Bromley et al. 2018; Du et al. 2018, 2019; Shen et al. 2018; Marchetti et al. 2019; Marchetti 2021; Marchetti et al. 2022; Liao et al. 2023). In addition to Gaia, large spectroscopic surveys have also proven to be powerful tools for discovering HVSs. For instance, Sun et al. (2025) recently identified ten new unbound B- and A-type HVS candidates using the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST).

The Sagittarius dwarf spheroidal galaxy (Sgr dSph) is the closest and third most massive satellite of the Milky Way, and it is currently undergoing significant tidal disruption (Ibata et al. 1994). The HVS and high-velocity star candidates originating from the Sgr dSph provide a unique window into its central dynamics and tidal history. Huang et al. (2021) reported the first such candidate, J1443+1453, a star with a galactocentric velocity of 599 km s⁻¹ that likely originated from the Sgr dSph. Building on this, Li et al. (2022) (hereafter Li22) combined proper motions from Gaia Early Data Release 3 (Gaia EDR3; Gaia Collaboration (2021)) with radial velocities from multiple spectroscopic surveys to identify 60 high-velocity candidates whose orbits are consistent with a Sgr origin. The locations of these stars on the Hertzsprung-Russell diagram, coupled with the chemical properties of the 19 candidates with measured metallic-ities, show similarities to members of the Sgr stream. The closest approach of the Li22 candidates to the Sgr dSph occurred during the progenitor’s pericentric passage (~38.2 Myr ago), leading the authors to suggest that these stars were tidally stripped from the system. Additionally, Li22 proposed that two HVSs within their sample might have originated from the Hills mechanism. Furthermore, Li et al. (2023) (hereafter Li23) assembled a large sample of 547 extreme-velocity stars with galactocentric velocities v_GSR > 0.8V_esc. Among these, 15 stars - including two HVS candidates - feature orbits that trace back to close encounters with the Sgr dSph. This hypothesized Sgr origin is further supported by an analysis of their distribution in the [α/Fe]-[Fe/H] plane.

Although many observational discoveries have been made, the dynamical mechanisms responsible for ejecting these stars have not been thoroughly explored. In this work, we investigate three possible origins for high-velocity stars linked to the Sgr dSph: the classic Hills mechanism, tidal stripping during the dwarf’s ongoing disruption, and chance passages of high-velocity halo stars through the Sgr dSph. The paper is structured as follows. In Section 2, we describe our selection of high-velocity stars from Gaia DR3, DESI DR1, and LAMOST DR12. In Section 3, we examine three possible origins for these stars: ejection via the Hills mechanism (Section 3.1), tidal stripping during the dwarf’s disruption (Section 3.2), and chance encounter of high-velocity halo stars with the Sgr dSph (Section 3.3). We provide a complementary chemical abundance analysis in Section 4. Finally, we present a concise discussion and our conclusions in Section 5.

2 Sample selection

2.1 Data

The Gaia DR3 catalog provides astrophysical parameters for approximately 470 million stars (Gaia Collaboration 2023) and radial velocities for 34 million sources in the Gaia DR3 Radial Velocity Spectrometer (RVS) measured sample (Katz et al. 2023). To ensure high-quality astrometric measurements, we applied the following criteria:

• (1)

  Renormalized unit weight error RUWE < 1.4. Small values mean the source is well fitted by a single star model.

• (2)

  ASTROMETRIC_EXCESS_NOISE_SIG < 2. Excess noise refers to the extra noise in each observation that causes the residual scatter in the astrometric solution. If ASTROMETRIC_EXCESS_NOISE_SIG is greater than two, then this excess noise cannot be ignored.

• (3)

  ASTROMETRIC_GOF_AL < 3. This parameter represents the goodness of fit between the astrometric model and the observed data. A higher value indicates a poor fit.

• (4)

  VISIBILITY_PERIODS_USED > 10. This parameter represents the number of visibility periods used in the astrometric solution. A high number of visibility periods is a better indicator of an astrometrically well-observed source.

• (5)

  ϖ - ϖ_ZP > 5 σ_ϖ. Here ϖ is the parallax and the σ terms denote its uncertainty. Parallax in the Gaia DR3 are known to exhibit a small systematic offset. We therefore corrected the reported parallaxes ϖ by subtracting the zero-point offset ϖ_ZP following the procedure of Lindegren et al. (2021).

• (6)

  For the Gaia DR3 RVS sample, we further included the following two criteria to ensure a robust radial velocity measurement: RV_EXPECTED_SIG_TO_NOISE ≥ 5 and RV_NB_TRANSITS ≥ 10. These cuts remove sources with low signal-to-noise radial velocity measurements and sources with an insufficient number of RVS transits.

Next, coordinates of selected Gaia DR3 RVS sample stars were transformed to the Galactocentric Cartesian frame with their 6D phase information by Astropy (Astropy Collaboration 2022). The heliocentric distance d_⊙ is equal to 1/(ω -ω_ZP), the Sun’s position (x_⊙,y_⊙,z_⊙) = (8.122, 0, 0.0208) kpc (GRAVITY Collaboration 2018; Bennett & Bovy 2019) was set, and the velocity (U_⊙, V_⊙, W_⊙) = (11.1, 245, 7.25)kms⁻¹ (Schönrich et al. 2010; McMillan 2017) was adopted. Then we imposed a velocity selection criterion of total velocity v_GC > 300 km s⁻¹. By combining these velocity thresholds with the previous quality cuts, we expected to retain only genuine highvelocity stars.

The DESI DR1 provides spectra for over 18 million unique targets observed in the Main Survey between May 2021 and June 2022 (Collaboration 2025). To enhance the scientific usability beyond the core spectroscopic outputs, the DESI Science Collaboration has produced a suite of value-added catalogs (VACs) accompanying DR1 (Koposov et al. 2025). In this study, we employed the Milky Way Survey (MWS) catalog, which contains detailed analyses of stellar spectra produced by the MWS Working Group, including stellar parameter determinations and quality assessments (Koposov et al. 2025). We combined radial velocities and quality flags from the MWS catalog with astrometric parameters and their uncertainties from Gaia DR3 through cross-matching.

To ensure a robust kinematic sample, we required RVS_WARN = 0, which indicates that the stars have reliable stellar model fits by the RVS pipeline. Finally, following the procedure applied to the Gaia DR3 RVS sample, we transformed the DESI DR1 MWS stars into Galactocentric coordinates and imposed a selection criterion of total velocity v_GC > 300 km s⁻¹.

LAMOST is a four-meter-aperture reflecting Schmidt telescope designed for efficiently obtaining spectra of both Galactic and extra-galactic objects in the optical band (Cui et al. 2012; Zhao et al. 2012). We utilized the A-, F-, G-, and K-type stellar parameter catalogs from the LAMOST DR12 low-resolution survey (LRS), which comprise a total of 8 370 041 spectra. Of these, 471 029 are A-type spectra, while the F-, G-, and K-type stars comprise 2 412 434; 3 966183; and 1520 395 spectra, respectively. Specifically, we used the radial velocities in LAMOST LRS DR12 derived by the LAMOST Stellar Parameter Pipeline (LASP; Luo et al. 2015), and we crossmatched other astrometric parameters with Gaia DR3. We applied the same quality filters as for the Gaia DR3 sample to the crossmatched LAMOST DR12-Gaia DR3 data and retained only sources with galactocentric total velocities υ_GC > 300 km s⁻¹.

2.2 Orbital analysis

To properly account for measurement uncertainties, we adopted a Bayesian framework. The posterior probability of the parameters θ = (d, v_α, v_δ, v_r) given the observables x = (ω, μ^*_α, μ_δ, rυ)^T is proportional to the product of the likelihood, which is modeled as a three-dimensional Gaussian distribution, and the prior probability: $$P(\theta \mid x)\propto \exp [-\frac{1}{2}(x-m(\theta ){)}^T{C}_x^{-1}(x-m(\theta ))]\times P_{\mathrm{p}\mathrm{r}\mathrm{i}\mathrm{o}\mathrm{r}}(d|\alpha ,\beta ,L),$$(1)

where $$m(\theta )=(1/d,v_{\alpha }/(kd),v_{\delta }/(kd),v_r{)}^T,k=4.74047,$$

and C_x is the covariance matrix of the observables. For the distance prior P_prior(d∖α,β, L), we used the geometric distance prior of Bailer-Jones et al. (2021), which is a three-parameter generalized gamma distribution and parameters are determined by each star’s HEALPix index.

We implemented Markov chain Monte Carlo sampling using the Python package emcee (Foreman-Mackey et al. 2013) to draw samples from the posterior distribution, generating 5000 Monte Carlo (MC) realizations for each star. Next, each realization’s coordinates and velocities were transformed to a Galactocen-tric frame. For each star, we took the median (50th percentile) of these realizations as the central estimate, with uncertainties defined by the 16th and 84th percentiles. In this work, we define high-velocity stars as those with galactocentric velocities greater than 400 km s⁻¹, and HVSs as those with galactocentric velocities exceeding 600 km s⁻¹. In previous studies, HVSs were commonly defined based on their escape probability (Marchetti et al. 2019; Marchetti 2021; Liao et al. 2023), determined by comparing a star’s velocity with the local escape velocity; stars with escape probabilities greater than 50% were classified as HVSs. We argue that this definition makes the classification strongly dependent on distance uncertainties, since stars with larger distance errors are more likely to be assigned high escape probabilities. Instead, we adopted a purely velocity-based definition, which we consider more robust. Subsequently, we compare our HVS candidates with the escape velocity curve to assess whether they are truly unbound.

For the high-velocity star candidates, we integrated their MC orbital realizations backward in time within a timedependent Galactic potential. This potential consists of the static MWPotential2014 (Bovy 2015) together with the moving Large Magellanic Cloud (LMC) and the moving Sgr dSph.

The LMC was modeled by a Hernquist potential (Hernquist 1990) with a total mass of M_LMC = 1.5 × 10¹¹ M_⊙ and a scale radius of r_LMC = 17.14 kpc (Erkal et al. 2019; Erkal & Belokurov 2020). At t = 0 (present day) its center of mass is placed at (α, δ) = (78.76°, −69.19°) with a heliocentric distance of 49.59kpc (Zivick et al. 2019; Pietrzynski et al. 2013), a proper motion of (μ*_α,μ_δ) = (1.91,0.229) mas, yr⁻¹, and a heliocentric line-of-sight velocity of 262.9 km s⁻¹ (van der Marel et al. 2002; Kallivayalil et al. 2013). We integrated the LMC orbit 1 Gyr into the past in MWPotential2014, including dynamical friction following Chandrasekhar (1943), and represented its time-varying contribution with the MovingObjectPotential built in galpy.

We modeled the Sgr dSph using a Plummer potential. To satisfy the observational constraints from Vasiliev & Belokurov (2020), who reported a mass of 4 × 10⁸ M_⊙ within 5 kpc and a peak circular velocity of 21 km s⁻¹ at ~2.5-3 kpc, we adopted the structural parameters from Rostami Shirazi et al. (2024): M_Sgr = 4.8 × 10⁸ M_⊙, a scale radius of r_sc = 1.8 kpc, and a half-mass radius of R_h = 2.34 kpc. The present-day Sgr dSph remnant is positioned at (α, δ) = (283.76°, −30.48°), with proper motions of (μ^*_α,μ_δ) = (−2.70, −1.35) mas yr⁻¹ and a heliocentric line-of-sight velocity of 142 km s⁻¹ (Vasiliev & Belokurov 2020; Vasiliev et al. 2021). The distance to the Sgr center remains the primary source of kinematic uncertainty, with values in the literature ranging from 24 to 28 kpc (Ibata et al. 1997; Hamanowicz et al. 2016; Ferguson & Strigari 2020; Vasiliev & Belokurov 2020; Vasiliev et al. 2021). Our tests indicate that while a smaller distance increases the number of high-velocity stars undergoing close encounters (as detailed in Section 3.3), the choice of distance does not significantly skew the resulting kinematic distributions (e.g., encounter velocities or flight times). Consequently, we adopted a distance of 27 kpc (Vasiliev et al. 2021), as it provides an optimal fit to the observed stream properties. We integrated the Sgr dSph orbit backward in time for 1 Gyr, accounting for dynamical friction and the gravitational influence of both the Milky Way and the LMC.

We quantified the probability that a given star experiences a close encounter with the Sgr dSph, denoted as P_Sgr, by $$P_{\mathrm{S}\mathrm{g}\mathrm{r}}=\frac{N_{\mathrm{M}\mathrm{C}}(r_{min}<2.5R_h)}{N_{\mathrm{M}\mathrm{C}}},$$(2)

where r_min is the closest distance between the star and the Sgr dSph along its backward-integrated orbit, N_MC = 5000 is the number of MC realizations used in this study, and N_MC(condition) is the number of MC realizations that satisfy the corresponding condition. The threshold 2.5R_h corresponds to two and a half times the remnant’s half-mass radius, providing a physically motivated boundary. In Gaia DR3, we identify 60 stars with P_Sgr > 0.5 out of 2967 high-velocity stars (2.02%), including six sources also reported by Li22. DESI DR1 contributes 17 such stars among 732 (2.32%), of which one overlaps with Li23. LAMOST DR12 contributes 18 such stars among 524 (3.44%), including five in common with Li23. These stars are displayed in Figure 1. We find only one star with v_GC exceeding 600 km s⁻¹ ; however, its galactocentric distance is about 2.25 kpc, and this velocity lies below the local escape velocity curve. We therefore conclude that no HVS originating from the Sgr dSph are identified in this work.

For each high-velocity star with P_Sgr > 0.5, we computed the relative velocity with respect to the Sgr dSph at pericentric passage (∆V) and the corresponding flight time to the present (T_f). The results are illustrated in Figure 2.

A clear clustering is observed for all three samples around T_f ~ 30-70 Myr and ∆V ~ 400-600 km s⁻¹, consistent with the flight time of the Sgr dSph (~39.7 Myr) since its most recent pericentric passage (Li et al. 2023). However, distinct differences appear in the extended distributions. Our sample exhibits a significant population extending towards long flight times (T_f ≳ 100 Myr) and correspondingly low relative velocities (∆V < 300 km s⁻¹). We interpret these outliers as halo stars, which is discussed further in Section 3.3. In comparison, Li22 shows a similar but sparser trend in this high-T_f/1ow-ΔV regime, while Li23 is tightly confined to the region associated with the recent Sgr passage.

Since the high-velocity stars in our sample experienced recent close encounters with the Sgr dSph, their motions are expected to be intrinsically linked to the dynamics of the Sgr system. To better characterize these kinematic patterns and facilitate a direct comparison with the Sgr stellar stream, we transformed the heliocentric coordinates of our candidates into the Sgr stream coordinate system (Λ_⊙, B_⊙). Given the variety of conflicting conventions in the literature, we adopted the right-handed heliocentric system as defined in Vasiliev et al. (2021) to ensure consistency with standard Galactic and International Celestial Reference System frames. Specifically, we followed the convention of Belokurov et al. (2014) for the longitudinal coordinate Λ_⊙, which represents the angle along the Sgr orbital plane. In this system, Λ_⊙ increases toward the leading arm (Λ_⊙ > 0) and decreases toward the trailing arm (Λ_⊙ < 0). For the latitude B_⊙, defined perpendicular to the orbital plane, we followed the convention of Majewski et al. (2003), which ensures that the orbital pole of the Sgr plane lies at B_⊙ = 90° (corresponding to (l, b) = (273.75°,-13.46°) in Galactic coordinates, or (α, δ) = (123.65°, −53.52°) in the International Celestial Reference System). In this coordinate system, the Sgr dSph remnant is centered at (Λ_⊙,0, B_⊙,0) ≈ (0°, 1.5°). We plot those high-velocity stars in Figure 3. In the left panel, the line-of-sight velocities v_r of these high-velocity stars display a clear correlation with the stream longitude Λ_⊙. The distribution is symmetric about Λ_⊙ = 50° : for Λ_⊙ > 50°, v_r > 0 (receding from the Sun), whereas for Λ_⊙ < 50°, v_r < 0 (approaching the Sun). This behavior arises from projection effects along the line of sight, since the stars share a similar orbital phase near pericentric passage with Sgr dSph - some having passed the Solar azimuth produce positive v_r, while others moving toward the Sun produce negative values.

In the right panel, the tangential velocity component $v_{{\mathrm{\Lambda }}_{\odot }}$ measured along the stream longitude exhibits a cosine-like variation with Λ_⊙, reaching its maximum at Λ_⊙ = 50° where the line-of-sight component crosses zero. This sinusoidal pattern likewise reflects geometrical projection of the common star motion onto the tangential direction. All three samples (this work, Li22, and Li23) reveal remarkably consistent trends, supporting the robustness of our selection procedure.

Fig. 1 Distribution of galactocentric total velocity (vGC) as a function of galactocentric radius (rGC) for high-velocity stars (PSgr > 0.5). Solid, dashed, and dotted curves show the escape-velocity profiles predicted by three Milky Way potential models: Bovy (2015), McMillan (2017), and Cautun et al. (2020).

Fig. 2 Relative velocity ∆V with respect to the Sgr dSph at pericentric passage versus flight time Tf on a logarithmic scale. Blue circles, orange triangles, and magenta squares represent results from this work, Li22, and Li23, respectively. Vertical red line indicates the 39.7 Myr flight time of the Sgr dSph since its most recent pericentric passage, with the shaded region representing the 1σ confidence interval.

3 Modeling possible origins of the Sgr high-velocity star candidates

3.1 The Hills mechanism

In this subsection, we investigate the hypothesis that HVSs can be produced at the center of the Sgr dSph via the Hills mechanism. We employed the Python package speedystar¹ (Evans et al. 2022b) to generate HVS samples from the Hills mechanism, releasing them at the corresponding orbital positions of the Sgr dSph, integrating their orbits within the Galactic potential and performing mock observation. Our study focuses on exploring a parameter space defined by three key variables: the assumed mass of the central black hole in the Sgr dSph (M_BH), the initial mass function (IMF) slope (κ), and the HVS ejection rate (η).

It is well established that the stellar velocity dispersion (σ) and the mass of a central black hole (M_BH) are strongly correlated (Kormendy & Ho 2013; Greene et al. 2020). Although the Sgr dSph is undergoing violent tidal disruption and lacks a classical bulge, its remnant velocity dispersion remains remarkably constant at σ ≈ 12-14 km s⁻¹ (Vasiliev & Belokurov 2020). To estimate the black hole mass, we considered two versions of the M_BH-σ relation: the canonical relation inferred from a large sample of galaxies (Kormendy & Ho 2013), where M_BH ∝ σ^4.4, and a version specifically tailored to low-mass black holes (Xiao et al. 2011), where M_BH ∝ σ^3.3. These relations yield a putative IMBH mass for the Sgr dSph within the range of 10³−10⁴ M_⊙. Furthermore, as shown by the analysis of the NewHorizon simulation (Dubois et al. 2021) in Beckmann et al. (2023), dwarf galaxies with stellar masses similar to that of the Sgr dSph are expected to possess IMBHs that remain in their seed state (M_BH < 10⁴ M_⊙). Consequently, we varied the parameter M_BH using values of [10³,5 × 10³,10⁴] M_⊙ to encompass the physically plausible range for such a system.

There have been several independent lines of evidence suggesting that the IMF in the Sgr dSph is top-light. From chemical abundance analysis of stars in the Sgr dSph, McWilliam et al. (2013) identified a ~0.4 dex deficiency in hydrostatic α-elements (O and Mg) relative to explosive α-elements (Si, Ca, and Ti), concluding that Sgr lacked enrichment from stars with masses ≳ 30 M_⊙. Furthermore, chemical evolution modeling of Sgr stream stars by Carlin et al. (2018) suggests that a steeper powerlaw index (κ ≈ 2.75) is necessary to match observed abundance patterns. Consequently, to encompass both the standard Salpeter slope and the steeper values suggested by these empirical constraints, we adopted a power-law form f(m) ∝ m^−κ with κ = [2.35,2.55,2.75,2.95,3.15,3.35].

The HVS ejection rate, η, is a relatively unconstrained parameter in our model. Theoretical estimates of the classical Hills mechanism in the GC typically range from 10⁻⁵ to 10⁻³ yr⁻¹ (Hills 1988; Yu & Tremaine 2003), depending on the rate at which stars are scattered into the “loss cone” - the region of phase space characterized by low-angular momentum orbits around Sgr A^*. Subsequent detailed simulations (e.g., Zhang et al. 2013), as well as calibrations based on the known population of HVS candidates (Bromley et al. 2012; Brown et al. 2014; Marchetti et al. 2018) and the frequency of tidal disruption events in the local Universe (Bromley et al. 2012; Brown 2015), yield consistent results. However, unlike the Milky Way center, the Sgr dSph is a system far from dynamical equilibrium due to ongoing tidal disruption. Since a self-consistent determination of the ejection rate through detailed dynamical simulations is beyond the scope of this work, we explored a wide range of values, η ∈ [10⁻², 10⁻³, 10⁻⁴, 10⁻⁵, 10⁻⁶] yr⁻¹, to encompass the plausible extremes and account for the significant uncertainties inherent to this disrupted system.

Other binary properties - such as the mass-ratio distribution f(q) ∝ q^γ and the semimajor axis distribution f (a) ∝ a^α - were shown by Evans et al. (2022a) to have negligible impact on the number of observable HVSs, and our experiments confirm insensitivity to γ and α in the Sgr dSph Hills scenario. Therefore, the default value γ = 0 and α = −1 are taken by recommendation from speedystar. We also must specify a metallicity distribution for the ejected stars - though our tests indicate negligible sensitivity of the detectable sample size to metallicity. Since there is currently no strong observational or theoretical evidence (Wrobel et al. 2011; Beckmann et al. 2023) indicating that the putative black hole at the center of the Sgr dSph must reside precisely at its dynamical center and coincide with the globular cluster M54, we adopted the metallicity distribution of the Sgr dSph core from Mucciarelli et al. (2017) for our simulation. This distribution was modeled as a normal distribution with <[Fe/H]) = −0.52 dex and σ = 0.17 dex.

Using a Monte Carlo approach with the above parameters, we randomly generated HVSs and assigned each a flight time uniformly distributed between 0 and 100 Myr. We excluded stars with flight time exceeding 100 Myr, as they travel beyond observable distances. We further accounted for prior main-sequence evolution by assuming each HVS has already consumed a portion of its main-sequence lifetime before ejection. Any star whose remaining main-sequence lifetime is shorter than its flight time was excluded from the simulation (we limited our analysis to unevolved stars and did not consider the evolved population for detection).

The above procedures yielded an HVS sample with known present-day positions, velocities, radii, luminosities, and effective temperatures. We computed their absolute magnitudes in the Gaia photometric bands and astrometric uncertainties following the method of Evans et al. (2022b), based on each star’s effective temperature, luminosity, and surface gravity, and accounting for line-of-sight extinction using the Galactic dust map combined15 of mwdust² (Bovy et al. 2016). We then estimated apparent magnitudes in the Gaia and Johnson-Cousins systems via the MIST model grids (Choi et al. 2016; Dotter 2016), and derived G_RVS magnitudes from polynomial fits of Jordi et al. (2010) using Gaia G and Johnson-Cousins V and I_c bands.

In Figure 4, we show the distribution of total velocities of Hills mechanism ejected HVS from a putative black hole at the center of the Sgr dSph, plotted in the Galactic-center rest frame as a function of galactocentric distance. In this simulation, we adopt a central black hole mass of 10⁴ M_⊙, an HVS ejection rate of 10⁻⁴yr⁻¹, and an initial mass function given by the classical Salpeter law (Salpeter 1955) with κ = 2.35. The resulting HVS span a very broad radial range - from the Solar neighborhood out to beyond 100 kpc - with most stars exceeding the local escape-velocity curve. This implies that, if HVSs are produced by the Hills mechanism in the Sgr dSph, their velocities would likely lie well above the escape speed, reaching on the order of 10³km s⁻¹. We also color-code each HVS by its Grvs magnitude, since this band places the most stringent observational limits: Gaia DR3 RVS sample requires G_RVS < 14, whereas Gaia DR4 RVS sample will extend to G_RVS < 16.2. Under this model, most simulated HVSs are too faint in the RVS band to be detected by Gaia DR3, and only a small portion would be accessible to Gaia DR4.

In Gaia DR3 and DR4, an HVS must satisfy G < 20.7 to be detectable. Since Gaia’s RVS sample remains the most efficient channel for identifying potential HVSs from the Sgr dSph (60 out of 95 candidates in this work are from the Gaia DR3 RVS sample), it is important to estimate how many HVSs can be included in the Gaia RVS sample. In the following, when we refer to an HVS as “observable” in Gaia DR3/DR4, we mean that it can be included in the corresponding RVS sample. Membership in the DR3 RVS sample further requires G_RVS < 14 for stars cooler than 6900 K or G_RVS < 12 for stars cooler than 14 500 K. In DR4, the requirements are relaxed to G_RVS < 16.2 for stars cooler than 6900 K or G_RVS < 14 for stars cooler than 14 500 K.

To estimate astrometric errors, we adopted the fivedimensional covariance distribution function for DR3 measurements proposed by Everall et al. (2021), and projected anticipated DR4 improvements by scaling down the DR3 errors. Radial-velocity uncertainties were estimated using PyGaia³, based on each star’s effective temperature and G_RVS. To distinguish HVSs from field stars, we further required each star’s total velocity to exceed the local escape velocity. Finally, to ensure precise orbit determinations, we demanded that the parallax uncertainty be less than 20% of the measured parallax.

For each parameter combination, we performed 100 simulations and applied 1000 bootstrap resamplings. For each bootstrap resample, we computed the mean number of observed HVSs, N_HVS. The median (50th percentile) of these resampled means was adopted as the central value, with uncertainties given by the 16th and 84th percentiles.

We find that N_HVS remains nearly zero across all parameter sets when considering Gaia DR3 observational constraints, as well as for cases with η ≤ 10⁻⁵ yr⁻¹ in the Gaia DR4 mock catalog. Consequently, Figure 5 only presents the results for the Gaia DR4 sample with η ≥ 10⁻⁴ yr⁻¹, where the predicted number of HVSs is statistically significant.

In DR4, N_HVS decreases significantly with increasing IMF slope κ. This trend reflects the transition from a canonical Salpeter to a top-light IMF, which reduces the number of highmass stars bright enough to be detected in DR4. We also find a modest increase in N_HVS with black hole mass, likely because a more massive black hole imparts higher ejection velocities, allowing HVSs to reach the solar neighborhood. Moreover, changing η by an order of magnitude shifts N_HVS by a similar factor. Our results agree with those of Evans et al. (2022a): only when η is high and the IMF is not top-light can a detectable sample be obtained. By adopting the proposed IMF slope of κ = 2.75 ± 0.2 (Carlin et al. 2018), we suggest that it is unlikely to detect genuine HVSs in Gaia DR4 produced via the Hills mechanism, except under conditions of extremely high ejection rates (η ≥ 10⁻²yr⁻¹).

Fig. 3 Distribution of high-velocity stars in Sgr Stream coordinates (Vasiliev et al. 2021). Left panel: radial velocity over Sgr stream coordinates Λ⊙. Right panel: Tangential velocity in the Λ⊙ direction over Λ⊙. Blue circles, orange triangles, and magenta squares represent results from this work, Li22, and Li23, respectively.

Fig. 4 Distribution of galactocentric total velocity (vGC) versus galac-tocentric radius (rGC) for HVSs ejected via the Hills mechanism from the center of the Sgr dSph, color coded by each star’s apparent magnitude in the Gaia RVS band. Results generated using the speedystar code assuming MBH = 104M⊙, ejection rate η = 10−4yr−1 and a Salpeter IMF (Salpeter 1955) with κ = 2.35. Solid, dashed, and dotted curves represent the escape-velocity profiles predicted by three Milky Way potential models (Bovy 2015; McMillan 2017; Cautun et al. 2020).

3.2 Tidal disruption of the Sgr dSph

In this subsection, we discuss whether the tidal disruption of the Sgr dSph could be a possible channel to produce observed highvelocity stars that have a close encounter with the Sgr dSph. In the pioneering work, Abadi et al. (2009) first proposed that stars stripped from a disrupted satellite galaxy could achieve velocities comparable to those produced by the Hills mechanism. Building upon this idea, Piffl et al. (2011) studied the scenario in detail by performing a suite of N-body simulations of a satellite’s encounter with a Milky Way-type host. They gathered statistics on the properties of stripped stars and showed that the maximum energy kick imparted to tidal debris depends primarily on the satellite’s initial mass. Their results indicate that only massive satellites (M_sat > 10⁹ M_⊙) can produce unbound stellar populations via tidal stripping.

To predict the kinematics of stripped stars, Piffl et al. (2011) developed a simple analytical model to estimate the maximum energy imparted to stars during the tidal disruption of the satellite. Applying this model to the Sgr dSph, which has a pericentric speed V_peri ≈ 350 kms⁻¹, an escape speed v_esc(R_tidal) ≈ 40 km s⁻¹, and an orbital angular momentum L_sat ≈ 5400 kpc km s⁻¹ (Vasiliev et al. 2021), one obtains ∆E_max ≈ 9000 km² s⁻². This value is much smaller than the magnitude of Sgr dSph’s orbital energy E_sat ≈ −4.9 × 10⁴ km² s⁻². To further validate the results presented above, we have constructed an analytical model that captures the detailed physical processes governing stellar stripping from the Sgr dSph. In this model, stars initially bound to Sgr dSph orbit its center with speed v_ini. As the satellite galaxy plunges through the Milky Way’s potential, its stars experience impulsive energy injections from the Galactic tidal field. Once a star acquires sufficient energy, it will escape from the Sgr dSph with a final speed of $${\mathbf{v}}_{\mathrm{f}\mathrm{i}\mathrm{n}}={\mathbf{v}}_{\mathrm{i}\mathrm{n}\mathrm{i}}+\mathrm{\Delta }\mathbf{v},$$(3)

where ∆v represents the velocity impulse imparted by the tidal interaction. For a fixed magnitude of ∆v, the magnitude of the final velocity |v_fin| is maximized if and only if ∆v is parallel to v_ini. To quantify ∆v, we computed the energy input using the so-called “tidal shock” formalism. In this framework, the change in a star’s kinetic energy is determined by integrating the timevarying tidal acceleration it experienced along the Sgr dSph’s orbit. We discuss the tidal shock approximation below.

The tidal shock formalism rests on two fundamental hypotheses. The first, known as the “impulse approximation”, posits that the internal velocity dispersion of the subject system is negligible compared to its relative velocity at pericenter passage. In the case of the Sgr dSph, the velocity dispersion, σ_v ≈ 13 km s⁻¹ (Vasiliev & Belokurov 2020), is indeed much smaller than its pericentric speed around the Milky Way, v_peri ≈ 350 km s⁻¹. Under this approximation, during pericenter passage (when tidal force is strongest) the stars within Sgr dSph can be treated as effectively stationary with respect to the galaxy’s center, allowing the tidal force to act impulsively.

The second fundamental hypothesis, known as the “distanttide approximation”, states that when the separation between the perturber (the Milky Way, in this case) and the subject system (the Sgr dSph) greatly exceeds the subject system’s characteristic size, the external potential varies smoothly across the system. Under this condition, one may expand the perturber’s gravitational field about the subject system’s center in a Taylor series. In a frame where the Sgr dSph’s center of mass is at the origin, let Φ(x, t) denote the gravitational potential of the MW. The tidal acceleration of star α relative to the origin is given by: $${\dot{\mathbf{v}}}_{\alpha }=-\sum _{j,k=1}^3{\hat{e}}_j{\mathrm{\Phi }}_{jk}{x}_{\alpha k},$$(4)

where $${\mathrm{\Phi }}_{jk}\equiv \frac{{\mathrm{\partial }}^2\mathrm{\Phi }}{\mathrm{\partial }{x}_j\mathrm{\partial }{x}_k}{|}_{\mathbf{x}=0}$$(5)

denotes the tidal tensor of the gravitational field induced by the perturber, evaluated at the origin of the Sgr-centric coordinate frame.

When both the distant-tide and high-speed (impulse) approximations hold simultaneously, the encounter is often referred to as a “tidal shock”. In the case of the Sgr dSph, the pericentric distance to the Milky Way is r_peri ≈ 16 kpc, while its tidal radius is r_t ≈ 2 kpc (Vasiliev & Belokurov 2020). Although the strict inequality r_t ≪ r_peri does not hold, we aim only to estimate the order of energy input from the tidal shock. Our estimated energy injection agrees well with the results of the N -body simulations presented by Vasiliev et al. (2021), lending support to the validity of this simplified approach. The physical meaning of Equation (4) is that the tidal acceleration on a star is correlated with the distance between it and the center of Sgr dSph, so we focus on stars at the tidal radius to evaluate the high-speed star generation by tidal disruption. For consistency with the preceding analysis, we adopted the MWPotential2014 model (Bovy 2015) to compute the tidal forces acting on the Sgr dSph. This potential is composed of a spherical bulge characterized by a power-law density profile with an exponential cut-off, an axisymmetric Miyamoto-Nagai disk, and a spherical Navarro-Frenk-White halo.

Suppose the trajectory of Milky Way in the Sgr-centric frame is given by the closed expression X(t). At any instant, one can perform a straightforward exercise to compute the tidal tensor at the center of mass of the Sgr dSph. The resulting tidal acceleration of a star α depends on time only through the trajectory X(t). We first computed the orbit of the Sgr dSph over the time interval [−0.5, 0.5] Gyr—encompassing its most recent pericentric passage around the Milky Way using galpy (Bovy 2015). Next, we performed spline interpolation on the orbit of Sgr dSph. We then numerically integrated the tidal acceleration of stars relative to the Sgr dSph centroid along this interpolated orbit to obtain the velocity increment ∆v. This computation was performed using the scipy library (Virtanen et al. 2020).

We randomly placed 1000 points at the tidal radius of the Sgr dSph and then integrated their velocity increments ∆v relative to the Sgr dSph’s center induced by the tidal forces. The results are shown in the left panel of Figure 6. We find that for most of the points at the tidal radius the velocity increment is approximately several tens km s⁻¹. According to Equation (3), if the direction of this increment aligns with the star’s initial velocity relative to the center of mass (typically 10-20 km s⁻¹; see Vasiliev & Belokurov 2020), stars attain the maximum escape velocity relative to the center, which is also below 150 km s⁻¹ . This is far smaller than the 400-600 km s⁻¹ relative velocities observed for the high-speed stars that have experienced a close encounter with the Sgr dSph.

These results can also be validated against the N-body model of Vasiliev et al. (2021), which provides a snapshot of the Sgr dSph following its tidal disruption by the Milky Way. We selected particles that were stripped within the last 100 Myr in this model and integrated their orbits backward to the moment they were tidally stripped. The outcome is presented in the right panel of Figure 6. We find that most of these particles approach within 2.5 half-mass radii of the Sgr dSph, justifying our use of this radius as the criterion for identifying stars that have undergone a close encounter. Moreover, all of these particles reach maximum velocities below 100 km s⁻¹ at pericenter, in good agreement with our estimates using the tidal-shock approximation. Therefore, for the Sgr dSph, tidal stripping alone is very unlikely to produce the high-velocity stars observed in our sample.

Fig. 5 Predicted number of observable HVSs, NHVS, in Gaia DR4 as a function of the IMF slope κ. Panels show results for different assumed central black hole masses in the Sgr dSph: MBH = 103M⊙ (left), 5 × 103M⊙ (center), and 104M⊙ (right). Solid orange, dashed sky-blue, and dotted bluish-green lines correspond to ejection rates of η = 10−2, 10−3, and 10−4yr−1, respectively. Shaded regions indicate the 16th-84th percentile range from bootstrap resampling.

Fig. 6 Left panel: distribution of the impulse imparted to stars uniformly distributed on the tidal radius of the Sgr dSph. Right panel: velocity difference between each star and the Sgr dSph’s center of mass at the stripping time, plotted against the normalized radial separation (∆r̄ = rmin/Rh), color coded by flight time.

3.3 Chance encounters of Milky Way halo stars with the Sgr dSph

In this subsection, we investigated the scenario in which highvelocity stars are simply Milky Way halo stars that naturally pass by the Sgr dSph at its pericenter before moving into the solar neighborhood. To assess this possibility, it is essential to construct a robust model of Milky Way background stars. For this purpose, we adopted the model developed by Binney & Vasiliev (2023) (hereafter BV23 model).

The adopted framework is defined by parameters governing the action-based distribution functions f(J) for four stellar disks (three thin-disk age cohorts and one thick disk), as well as for a spheroidal bulge and spheroidal components of both the stellar and dark halos. Given these distribution functions and a specified gas distribution, one can self-consistently solve for the spatial densities of stars and dark matter and the gravitational potential they generate. The first constraint of the BV23 model is the kinematics of stars from Gaia DR2 (Gaia Collaboration 2018a) that possess the Radial Velocity Spectrometer measured line-of-sight velocities (Gaia Collaboration 2018b). The stellar locations and velocities of RVS stars used in the BV23 model are computed by Schönrich et al. (2019). The second constraint is the vertical stellar density profile in the column above the Sun (Gilmore & Reid 1983). The model is implemented using the Agama code (Vasiliev 2019).

We first sampled 100 million stars from the BV23 model, then selected stars with heliocentric distances less than 5 kpc, as parallax measurements are more reliable within this range. Next, we applied the same velocity cut as in Section 2, identifying stars whose galactocentric total velocity exceeds 400 km s⁻¹, which left us with 2408 stars. We then selected the subset of these stars with r_min < 2.5R_h, leaving 81 stars in total. The fraction of these stars whose orbits trace back to the Sgr dSph relative to the total high-velocity sample is approximately 3.36%, which closely matches the observed fraction (2.02% for Gaia DR3 RVS, 2.32% for DESI DR1, and 3.44% for LAMOST DR12).

Figure 7 presents the galactocentric total velocity (v_GC) as a function of galactocentric radius (r_GC) for the simulated stars that satisfy the aforementioned selection criteria. The simulated sample traces the same locus as the bulk of the observed stars in Figure 1. Notably, none of the simulated stars exceeds the local escape speed predicted by any of the Galactic potentials adopted in this study (Bovy 2015; McMillan 2017; Cautun et al. 2020), which is because actions can only be computed for bound stars.

In Figure 8, we present the relative velocities at pericentric passage and the corresponding flight times for high-velocity stars in the BV23 model that undergo close encounters with the Sgr dSph. The left panel shows the distribution of relative velocities, while the right panel shows the distribution of flight times, and for comparison we overplot the observed distributions of high-velocity stars that encounter the Sgr dSph. The ranges of both relative velocity and flight time in the BV23 model closely match the observations, and the model also reproduces the “low-velocity tail” and the “long-flight-time tail.” This agreement indicates that the model captures the underlying physical processes in the data rather than merely reflecting selection effects.

In Figure 9, we show the distribution of high-velocity stars in the BV23 model that undergo close pericentric encounters with the Sgr dSph, expressed in the Sgr-stream coordinate system. The left panel plots radial velocity v_r as a function of Λ_⊙, while the right panel shows the tangential velocity along the stream longitude, v_Λ⊙, versus Λ_⊙. The BV23 model captures the main kinematic features seen in the observations, including the positive correlation between v_r and Λ_⊙ and the cosine-like variation of v_Λ⊙ with Λ_⊙.

To further quantify the statistical consistency between the observations and the BV23 model, we performed a series of nonparametric tests. Two-sample Kolmogorov-Smirnov (KS) tests yield p-values of 0.37 for the relative velocities and 0.36 for the flight times. Furthermore, a two-sample energy test (Székely & Rizzo 2013) yields p-values of 0.33 for the v_r-Λ_⊙ distribution and 0.10 for the v_Λ⊙-Λ_⊙ projection. These high p-values indicate that the observed and modeled distributions are statistically indistinguishable.

Combining the discussion in this subsection with the previous two, we find that for the Sgr dSph, the two mechanisms traditionally considered as potential channels for producing high-velocity stars - the Hills mechanism and tidal disruption - both face difficulties in explaining the observed stars whose orbits trace back to the Sgr dSph. High-velocity stars produced by the Hills mechanism are unlikely to be detected with current observational facilities, while tidal disruption is theoretically inefficient in generating such stars. By contrast, our analysis of the BV23 model shows that chance encounters of Milky Way halo stars with the Sgr dSph yield a distribution that closely resembles the observations. This suggests that the observed high-velocity stars are most likely halo stars of the Milky Way. Furthermore, this provides a natural explanation for why adopting a smaller distance to the Sgr center results in a higher count of candidate stars, as noted in Section 2.2. A smaller assumed distance implies that the Sgr dSph penetrates deeper into the Milky Way (reaching smaller galactocentric radii) during its pericentric passage. This deeper penetration into the denser regions of the Galactic halo facilitates a larger number of chance encounters with halo stars, thereby increasing the number of identified candidates.

It must be emphasized, however, that our discussion here is preliminary and strongly model-dependent. First, the BV23 model is an N-body simulation, and we have not performed mock observations to account for uncertainties in the phasespace coordinates of stars. Such a procedure requires consideration of the stellar mass distribution, metallicity, dust extinction, and telescope selection functions across the Galaxy, which is a highly nontrivial task that we plan to address in future work. Second, while the BV23 model provides a relatively robust description of the Galactic disk, it may be insufficient for the halo, where the influence of the Gaia-Sausage-Enceladus merger likely needs to be considered. Using other Galactic models could lead to different conclusions. For example, when we employed the widely used GEDR3 mock model (Rybizki et al. 2020), we found only about ten stars in the entire dataset whose orbits intersect the Sgr dSph. A detailed comparison and discussion of different Galactic models will be presented in future work.

Fig. 7 Distribution of galactocentric total velocity (vGC) versus galactocentric radius (rGC) for stars sampled from Binney & Vasiliev (2023) model with vGC > 400km s−1 and rmin < 2.5Rh. Solid, dashed, and dotted curves show the escape-velocity profiles predicted by three Milky Way potential models: Bovy (2015), McMillan (2017), and Cautun et al. (2020).

Fig. 8 Left panel: distribution of relative velocities at pericentric passage (∆V) for high-velocity stars encountering the Sgr dSph, with observations shown as solid blue lines and the BV23 model as dashed black lines. Right panel: distribution of flight times since pericenter (Tf). The line styles are the same as in the left panel.

Fig. 9 Distribution of high-velocity stars in the Sgr-stream coordinate system. Left panel: radial velocity vr as a function of Λ⊙. Right panel: tangential velocity along the stream longitude, vΛ⊙, versus Λ⊙. The black circles represent stars from the BV23 model, and blue dots indicate the observed high-velocity stars that encounter the Sgr dSph in this work.

4 Chemical analysis

In this section, we conducted a complementary chemical abundance analysis to further substantiate our findings. Among our 95 high-velocity star candidates, 17 from DESI DR1 and 18 from LAMOST DR12 have metallicity measurements. However, the [a/Fe] measurements from DESI exhibit data quality issues, as crossmatching with APOGEE DR17 reveals significant scatter and systematic biases. Furthermore, only a limited subset of the LAMOST sample has [a/Fe] measurements. Consequently, we restricted our chemical analysis to [Fe/H].

Figure 10 displays the metallicity distribution of our highvelocity candidates, which exhibits a bimodal structure with peaks at [Fe/H]~ −1.3 and approximately −2.1 dex. For comparison, we plot the metallicity distributions of the Sgr dSph and its stream from Hayes et al. (2020), which include 166 stream members and 710 members of the Sgr dSph core. The Sgr stream members peak at [Fe/H]~ −1, while the Sgr core members peak at approximately −0.5, reflecting the metallicity gradient of the progenitor prior to tidal disruption. We also compare our sample with the Milky Way halo population characterized by Liu et al. (2018) using 4680 giant stars (|z| > 5 kpc) from LAMOST. Their study describes the halo metallicity profile as a bimodal Gaussian with peaks at [Fe/H]~ −1.2 ± 0.3 (inner halo) and −2.0 ± 0.2 (outer halo), comprising 91 and 9% of the population, respectively.

The chemical composition provides a crucial test for the origin of these stars. Notably, the bulk of the Sgr stream was stripped during the most recent pericentric passage approximately 1 Gyr ago. Li22 and Li23 invoked the chemical similarity between high-velocity stars and the Sgr stream as evidence to support a Sgr origin. However, our high-velocity candidates experienced close encounters with the Sgr dSph within the last 100 Myr, a much shorter timescale than the typical stripping history of the stream. If these stars indeed originated from the Sgr dSph, they would have had to be ejected by a violent dynamical process within the core. Nevertheless, the metallicity of our high-velocity sample deviates significantly from that of the Sgr core.

Furthermore, a two-sample KS test yields a p-value of 5.38 × 10⁻⁷ when comparing the Sgr stream members to our sample, strongly rejecting the hypothesis that these two populations are drawn from the same distribution. As shown in Figure 10, while there is overlap in the metallicity range, the distribution peaks are distinct: the Sgr stream is metal-richer by ~0.2-0.3 dex compared to both our candidates and the Milky Way halo. Conversely, a two-sample KS test between our sample and the halo distribution yields a p-value of 0.019. Given that significance thresholds are statistical conventions without intrinsic physical meaning, the key point is that this p-value is orders of magnitude higher than the result for the Sgr stream (for comparison, the KS test between the Sgr stream and the Milky Way halo yields p = 7.35 × 10⁻¹⁴).

However, we identify three very metal-poor (VMP) stars with [Fe/H] below −2.5 dex that cannot be adequately represented either by the Sgr or Milky Way halo distributions. These stars have metallicities of −3.48 ± 0.08 dex, −3.99 ± 0.04 dex, and −2.71 ± 0.09 dex, warranting further discussion. As noted by Liu et al. (2018), their halo sample is relatively small and incomplete at the very metal-poor end. After excluding these three stars from our high-velocity sample, the p-value between the truncated sample and the Milky Way halo increases to 0.056, while the p-value between the truncated sample and the Sgr stream increases to 2.77 × 10⁻⁶.

All three VMP stars originate from DESI DR1, with G-band magnitudes of 19.88 mag, 19.15 mag, and 17.46 mag, respectively. Given that LAMOST observations are challenging at such faint magnitudes, these stars are not well represented by the Milky Way metallicity distribution from Liu et al. (2018), which is based on LAMOST data. In contrast, DESI provides an order-of-magnitude increase in the number of faint stars with 17.5 < G < 21 that have radial velocity and abundance measurements compared to existing spectroscopic surveys, making it highly suitable for identifying faint VMP stars (Koposov et al. 2025). While a more comprehensive Milky Way metallicity distribution extending to the very-metal-poor end could be constructed from DESI data, this investigation is beyond the scope of the present work.

Fig. 10 Metallicity ([Fe/H]) distribution of the 95 high-velocity star candidates compared with Sgr and Milky Way halo populations. The gray histogram shows our candidates. Sgr stream members (dark cyan shading) and Sgr dSph core members (purple shading) are from Hayes et al. (2020). Milky Way halo giant stars from Liu et al. (2018) are shown with a solid blue line, encompassing the inner (dashed green) and outer (dashed red) halo components.

5 Discussion and conclusions

In this work, we used Gaia DR3, DESI DR1, and LAMOST DR12 to identify 95 high-velocity stars with close encounters to the Sgr dSph. Their relative velocities (∆V ~ 400-600 km s⁻¹ ) and short flight times (T_f < 100 Myr) coincide with the epoch of Sgr’s recent pericentric passage, and their kinematics follow coherent trends in Sgr stream coordinates.

We systematically tested the three possible origins proposed in previous studies. Simulations show that the Hills mechanism would not yield detectable HVSs in Gaia DR3 and only marginally in DR4 under extreme assumptions, while tidal disruption provides velocity kicks (∆V < 150 km s⁻¹) far below the observed range. By contrast, the Binney & Vasiliev (2023) halo model reproduces the observed velocity and orbital distributions, indicating that most candidates are ordinary halo stars. This conclusion is further supported by our metallicity distribution analysis, which indicates that these high-velocity stars are chemically consistent with the Milky Way halo population.

There are other mechanisms capable of producing highvelocity stars that were not explicitly discussed in previous sections. First, a supernova origin can be confidently excluded because the gas in the Sgr dSph was stripped approximately 1 Gyr ago (Tepper-Garcia & Bland-Hawthorn 2018), consistent with the last substantial burst of star formation in Sgr. This time interval significantly exceeds the lifespan of massive stars that undergo core-collapse. Second, for the same reason, the Sgr dSph possesses no young massive clusters. Third, the Sgr dSph currently hosts four bound globular clusters: M54, Terzan 7, Terzan 8, and Arp 2, with M54 located coincident with the core of the Sgr dSph. While single star-binary interactions involving compact objects in globular clusters can eject high-velocity stars (Weatherford et al. 2023; Evans et al. 2025), the ejection rate for stars with velocities exceeding 400 km s⁻¹ (as observed in this work) is two or three orders of magnitude lower than that of the Hills mechanism (Huang et al. 2025; Evans et al. 2025), and is therefore negligible.

Thus, our results suggest that both the candidates identified in this work and those reported in earlier studies through orbital back-tracing to the Sgr dSph are most likely ordinary halo stars. This highlights the need for caution when inferring the origin of high-velocity stars solely from orbital reconstruction, as the contribution of halo stars must be carefully accounted for. Future studies may benefit from incorporating more detailed chemical abundance information to obtain more robust conclusions. Since halo stars are, by definition, bound to the Galaxy, a secure detection of truly unbound HVSs originating from the Sgr dSph would immediately point to intriguing and extreme dynamical processes - though no such sources are found in this study. While our analysis indicates that genuine Sgr-origin HVSs are unlikely to be identified under current observational limits, this work provides a useful stepping stone for future searches.

Acknowledgements

This work was supported by National Key R&D Program of China No. 2024YFA1611900, and the National Natural Science Foundation of China (NSFC Nos. 11973042, 11973052). 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 research used data obtained with the Dark Energy Spectroscopic Instrument (DESI). DESI construction and operations is managed by the Lawrence Berkeley National Laboratory. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of High-Energy Physics, under Contract No. DE-AC02-05CH11231, and by the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility under the same contract. Additional support for DESI was provided by the U.S. National Science Foundation (NSF), Division of Astronomical Sciences under Contract No. AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research Laboratory; the Science and Technology Facilities Council of the United Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation; the French Alternative Energies and Atomic Energy Commission (CEA); the National Council of Humanities, Science and Technology of Mexico (CONAH-CYT); the Ministry of Science and Innovation of Spain (MICINN), and by the DESI Member Institutions: www.desi.lbl.gov/collaborating-institutions. The DESI collaboration is honored to be permitted to conduct scientific research on I’oligam Du’ag (Kitt Peak), a mountain with particular significance to the Tohono O’odham Nation. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the U.S. National Science Foundation, the U.S. Department of Energy, or any of the listed funding agencies. LAMOST is a National Major Scientific Project built by the Chinese Academy of Sciences, which has been provided by the National Development and Reform Commission.

References

All Figures

	Fig. 1 Distribution of galactocentric total velocity (vGC) as a function of galactocentric radius (rGC) for high-velocity stars (PSgr > 0.5). Solid, dashed, and dotted curves show the escape-velocity profiles predicted by three Milky Way potential models: Bovy (2015), McMillan (2017), and Cautun et al. (2020).
In the text

Fig. 2 Relative velocity ∆V with respect to the Sgr dSph at pericentric passage versus flight time Tf on a logarithmic scale. Blue circles, orange triangles, and magenta squares represent results from this work, Li22, and Li23, respectively. Vertical red line indicates the 39.7 Myr flight time of the Sgr dSph since its most recent pericentric passage, with the shaded region representing the 1σ confidence interval.

In the text

	Fig. 3 Distribution of high-velocity stars in Sgr Stream coordinates (Vasiliev et al. 2021). Left panel: radial velocity over Sgr stream coordinates Λ⊙. Right panel: Tangential velocity in the Λ⊙ direction over Λ⊙. Blue circles, orange triangles, and magenta squares represent results from this work, Li22, and Li23, respectively.
In the text

Fig. 4 Distribution of galactocentric total velocity (vGC) versus galac-tocentric radius (rGC) for HVSs ejected via the Hills mechanism from the center of the Sgr dSph, color coded by each star’s apparent magnitude in the Gaia RVS band. Results generated using the speedystar code assuming MBH = 104M⊙, ejection rate η = 10−4yr−1 and a Salpeter IMF (Salpeter 1955) with κ = 2.35. Solid, dashed, and dotted curves represent the escape-velocity profiles predicted by three Milky Way potential models (Bovy 2015; McMillan 2017; Cautun et al. 2020).

In the text

Fig. 5 Predicted number of observable HVSs, NHVS, in Gaia DR4 as a function of the IMF slope κ. Panels show results for different assumed central black hole masses in the Sgr dSph: MBH = 103M⊙ (left), 5 × 103M⊙ (center), and 104M⊙ (right). Solid orange, dashed sky-blue, and dotted bluish-green lines correspond to ejection rates of η = 10−2, 10−3, and 10−4yr−1, respectively. Shaded regions indicate the 16th-84th percentile range from bootstrap resampling.

In the text

	Fig. 6 Left panel: distribution of the impulse imparted to stars uniformly distributed on the tidal radius of the Sgr dSph. Right panel: velocity difference between each star and the Sgr dSph’s center of mass at the stripping time, plotted against the normalized radial separation (∆r̄ = rmin/Rh), color coded by flight time.
In the text

	Fig. 7 Distribution of galactocentric total velocity (vGC) versus galactocentric radius (rGC) for stars sampled from Binney & Vasiliev (2023) model with vGC > 400km s−1 and rmin < 2.5Rh. Solid, dashed, and dotted curves show the escape-velocity profiles predicted by three Milky Way potential models: Bovy (2015), McMillan (2017), and Cautun et al. (2020).
In the text

	Fig. 8 Left panel: distribution of relative velocities at pericentric passage (∆V) for high-velocity stars encountering the Sgr dSph, with observations shown as solid blue lines and the BV23 model as dashed black lines. Right panel: distribution of flight times since pericenter (Tf). The line styles are the same as in the left panel.
In the text

	Fig. 9 Distribution of high-velocity stars in the Sgr-stream coordinate system. Left panel: radial velocity vr as a function of Λ⊙. Right panel: tangential velocity along the stream longitude, vΛ⊙, versus Λ⊙. The black circles represent stars from the BV23 model, and blue dots indicate the observed high-velocity stars that encounter the Sgr dSph in this work.
In the text

Fig. 10 Metallicity ([Fe/H]) distribution of the 95 high-velocity star candidates compared with Sgr and Milky Way halo populations. The gray histogram shows our candidates. Sgr stream members (dark cyan shading) and Sgr dSph core members (purple shading) are from Hayes et al. (2020). Milky Way halo giant stars from Liu et al. (2018) are shown with a solid blue line, encompassing the inner (dashed green) and outer (dashed red) halo components.

In the text