# A Pristine-UNIONS view on the Galaxy: Kinematics of the distant spur feature of the Sagittarius stream traced by blue horizontal branch stars

M. Bayer, E. Starkenburg, G F Thomas, N F Martin, A. Helmi, A. Byström, T. de Boer, E. Fernández Alvar, S. Gwyn, R. Ibata, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

M. Bayer, E. Starkenburg, G F Thomas, N F Martin, A. Helmi, et al., A Pristine-UNIONS view on the Galaxy: Kinematics of the distant spur feature of the Sagittarius stream traced by blue horizontal branch stars. Astronomy & Astrophysics - A&A, 2025, 701, pp.A117. 10.1051/0004-6361/202554131. insu-05263438

HAL Id: insu-05263438

https://insu.hal.science/insu-05263438v1

Submitted on 16 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Pristine-UNIONS view on the Galaxy: Kinematics of the distant spur feature of the Sagittarius stream traced by blue horizontal branch stars

M. Bayer<fnref n="1" />,*, D. E. Starkenburg<fnref n="1" />, G. F. Thomas<fnref n="2" />,<fnref n="3" />, N. F. Martin<fnref n="4" />,<fnref n="5" />, A. Helmi<fnref n="1" />, A. Byström<fnref n="6" />, T. de Boer<fnref n="7" />, E. Fernández Alvar<fnref n="2" />,<fnref n="3" />, S. Gwyn<fnref n="8" />, R. Ibata<fnref n="4" />, P. Jablonka<fnref n="9" />, G. Kordopatis<fnref n="10" />, T. Matsuno<fnref n="11" />, A. W. McConnachie<fnref n="8" />, G. E. Medina<fnref n="12" />,<fnref n="13" />, S. Rusterucci<fnref n="4" />,<fnref n="1" />, R. Sánchez-Janssen<fnref n="14" />, F. Sestito<fnref n="15" />, and A. Viswanathan<fnref n="16" />

Received 13 February 2025 / Accepted 8 July 2025

## ABSTRACT

**Context.** Providing a detailed picture of the Sagittarius stream offers important constraints on the build-up of the Galactic halo as well as its gravitational potential at large radii. While several attempts have been made to model the structure of the Sagittarius stream, no model has yet been able to match all the features observed for the stream. Moreover, for several of these features, the observational characterisation of their properties is rather limited, particularly at large distances.

**Aims.** The aim of this work is to investigate the kinematics of the Sagittarius stream’s outermost spur feature using blue horizontal branch (BHB) stars.

**Methods.** Candidate BHB stars were selected by combining two approaches: one capitalising on Pan-STARRS1 3P griz and u photometry taken as part of UNIONS, the other using Pristine Survey CaHK and SDSS ugr photometry. Follow-up optical spectra were obtained using ESO/VLT/FORS2 to confirm their BHB nature and obtain line-of-sight (LOS) velocities.

**Results.** Of our 25 candidates, 20 stars can be confirmed as bona fide BHB stars. Their LOS velocities, together with the 3D positions of these stars qualitatively match well with Sagittarius model predictions and trace the outer apocentre of the trailing arm and its spur feature very clearly. The quantitative offsets that are found between our data and the different models can be used to provide information about the Galactic gravitational potential at large distances. We present a first tentative analysis in this direction, and show that a realistic state-of-the-art N-body simulation of Sagittarius would provide better agreement with our observations if the enclosed mass of the Milky Way within 100 kpc were lowered to $(5.3 \pm 0.4) \times 10^{11} \text{ M}_\odot$ (versus $(5.6 \pm 0.4) \times 10^{11} \text{ M}_\odot$).

**Conclusions.** Our selection of bona fide BHB stars provides a new view on the outermost structure in 3D positions and LOS velocities of the Sagittarius debris.

**Key words.** techniques: radial velocities – techniques: spectroscopic – stars: horizontal-branch – stars: kinematics and dynamics – Galaxy: halo

## 1. Introduction

The Sagittarius dwarf spheroidal galaxy is one of the largest building blocks of the Milky Way (MW) halo (e.g. Majewski et al. 2003; Read & Erkal 2019; Naidu et al. 2020) with an

estimated present total mass inside the remnant of at least $10^9 \text{ M}_\odot$ (Ibata et al. 1997). Over the past 30 years since its original discovery (Ibata et al. 1994), it has become clear that the stellar stream produced by the infall of this dwarf galaxy onto the MW is wrapping around the Galaxy fully, and that it is a major component of the Galactic halo (e.g. Ibata et al. 2001;

<footnote n="1">¹ Kapteyn Astronomical Institute, University of Groningen, Landleven 12, 9747AD Groningen, The Netherlands</footnote>

<footnote n="2">² Instituto de Astrofísica de Canarias, 38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="3">³ Universidad de La Laguna, Dpto. Astrofísica, 38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="4">⁴ Université de Strasbourg, Observatoire astronomique de Strasbourg, UMR 7550, 67000 Strasbourg, France</footnote>

<footnote n="5">⁵ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="6">⁶ Institute for Astronomy, University of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="7">⁷ Institute for Astronomy, University of Hawaii, 2680 Woodlawn Drive, Honolulu, HI 96822, USA</footnote>

<footnote n="8">⁸ NRC Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada</footnote>

<footnote n="9">⁹ Laboratoire d'astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), 1290 Sauverny, Switzerland</footnote>

<footnote n="10">¹⁰ Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, 06000 Nice, France</footnote>

<footnote n="11">¹¹ Astronomisches Rechen-Institut, Zentrum für Astronomie der Universität Heidelberg, Mönchhofstraße 12–14, 69120 Heidelberg, Germany</footnote>

<footnote n="12">¹² David A. Dunlap Department of Astronomy & Astrophysics, University of Toronto, 50 St George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="13">¹³ Dunlap Institute for Astronomy & Astrophysics, University of Toronto, 50 St George Street, Toronto, ON M5S 3H4, Canada</footnote>

<footnote n="14">¹⁴ STFC UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK</footnote>

<footnote n="15">¹⁵ Centre for Astrophysics Research, Department of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK</footnote>

<footnote n="16">¹⁶ Department of Physics and Astronomy, University of Victoria, PO Box 3055, STN CSC, Victoria, BC V8W 3P6, Canada</footnote>

<footnote>* Corresponding author: mbayer@astro.rug.nl</footnote>
Majewski et al. 2003; Helmi 2004; Johnston et al. 2005; Belokurov et al. 2006, 2014; Sesar et al. 2017a; Antoja et al. 2020; Ibata et al. 2020; Ramos et al. 2020; Naidu et al. 2020; Ramos et al. 2022). Specifically, the stellar stream is believed to dominate the outer Galactic halo at galactocentric radii beyond 25 kpc (Naidu et al. 2020).

Many intricate features can be seen in the arms of the Sagittarius debris. With the onset of large-area deep photometric surveys such as the Sloan Digital Sky Survey (SDSS, York et al. 2000), positional and photometric data of giant stars with relatively precise distances, including stars on the horizontal branch (HB), were used to map the streams. It was established that the apocentres of the leading and trailing arm are significantly different (47.8 ± 0.5 and 102.5 ± 2.5 kpc, respectively, Belokurov et al. 2014). Belokurov et al. (2014) furthermore reported (i) a confirmation of the bifurcation of the southern trailing tail (Koposov et al. 2012) from the double-peaked profile of blue HB (BHB) stars across the stream at longitude 145° < Λ<sub>Sgr</sub> < 170° in stream coordinates<fnref n="1" />, (ii) a gentle precession of the orbital plane, and (iii) a potentially differing evolution across the sky of the brighter and fainter tails.

Further insight in the 3D structure of the Sagittarius debris, out to large distances, can be provided by variable RR Lyrae stars (e.g. Vivas et al. 2005). These tracers can be uncovered in time-domain photometry, as for instance provided by the Panoramic Survey Telescope and Rapid Response System 1 (Pan-STARRS 1, hereafter PS1, Chambers et al. 2016), or currently with much larger coverage and precision by Gaia DR3 (Gaia Collaboration 2023; Clementini et al. 2023; Li et al. 2023b; Muraveva et al. 2025). The PS1 dataset resulted in the remarkably clean and complete catalogue of RR Lyrae stars by Sesar et al. (2017b) and Hernitschek et al. (2017). Sesar et al. (2017a) use this sample within 13° of the Sagittarius orbital plane to study the Sagittarius debris, resulting in a selection of 19 000 candidate RR Lyrae stars. With their catalogue sensitivity out to larger distances, they reveal an extra arm at heliocentric distances of 120 kpc in the trailing stream at Λ<sub>Sgr</sub> = 172° called the spur feature. Starkenburg et al. (2019) subsequently note that their sample of photometrically selected candidate BHB stars also shows the same feature, confirming that the most distant features of the Sagittarius stream are equally complex, with many intricate features.

Many studies have demonstrated they can model the multiple
wraps of the Sagittarius debris (e.g. Ibata & Lewis 1998; Helmi
& White 2001; Law & Majewski 2010; Dierickx & Loeb 2017;
Thomas et al. 2017; Fardal et al. 2019; Vasiliev et al. 2021; Oria
et al. 2022; Davies et al. 2024; Hainje et al. 2025). However, only
a few studies have made predictions for the spur feature under
study. For instance, both Dierickx & Loeb (2017) and Fardal et al.
(2019) show that by using N-body simulations one can repro-
duce the multiple wraps as well as the spur feature. As the data
are consistent with a wide range of models, the authors mention
several tests for their models, including the velocity structure of
the spur feature. Fardal et al. (2019) provide a number of impor-
tant model predictions to separate stars dynamically in the spur
feature from the main trailing arm. These authors predict that
line-of-sight (LOS) velocity measurements for a sample of bright
spur feature members reaching an uncertainty of about 10 km s⁻<fnref n="1" />
would be sufficient for this purpose.

An additional study to predict the stream feature, and fur-
thermore also model the influence of the LMC on the merger,

as first suggested by Vera-Ciro & Helmi (2013), is presented by
Vasiliev et al. (2021). Their N-body simulations of the interac-
tion of the LMC, our Galaxy, and Sagittarius produce a spur
feature composed of stars stripped about 2.5 billion years ago
at the penultimate pericentre passage of the Sagittarius galaxy.
Subsequently, Oria et al. (2022) re-explore the idea that the pro-
genitor of Sagittarius was a disc galaxy, with the disc roughly
perpendicular to both the MW disc and the Sagittarius dwarf
galaxy orbital plane. Their model results in the observed bifur-
cation of both the leading and trailing part of the Sagittarius tidal
arm (as observed in more detail by, e.g. Ramos et al. 2022), but it
does not seem to significantly affect the Sagittarius spur region.

In a recent study, Hainje et al. (2025) proposed a novel
hypothesis concerning the correlation of dark matter (DM) self-
interactions and the observed spur feature. The authors note that
the presence of DM with self-interactions of sufficient magni-
tude can effectively remove mass from the progenitor of the
Sagittarius dwarf galaxy via the scattering of DM particles in
the Sagittarius dwarf galaxy and the Galaxy (a phenomenon
they term ram pressure evaporation). In this way, self-interaction
can result in a significant increase in the loss of stellar and DM
mass from the progenitor (see also their figure 2). Consequently,
the enhanced mass loss of the progenitor due to ram pressure
evaporation can dramatically suppress the spur, showcasing how
observational constraints on the spur feature might ultimately
present constraints on dark matter properties.

In addition to the more diffuse stellar wraps of the Sagittarius streams, several globular clusters are associated with its disruption, including NGC 2419, which is close in the sky and in distance to the spur feature (e.g. Irwin 1999; Newberg et al. 2003; Ruhland et al. 2011; Belokurov et al. 2014; Sohn et al. 2018; Massari et al. 2019; Antoja et al. 2020; Bellazzini et al. 2020; Peñarrubia & Petersen 2021; Davies et al. 2024; Chen & Gnedin 2024; Rostami Shirazi et al. 2024); for example, previous studies have used the phase space (e.g. Irwin 1999; Newberg et al. 2003; Ruhland et al. 2011; Belokurov et al. 2014; Sohn et al. 2018; Antoja et al. 2020; Peñarrubia & Petersen 2021; Davies et al. 2024; Rostami Shirazi et al. 2024) plus age constraints (e.g. Massari et al. 2019) or chemical abundances (e.g. Bellazzini et al. 2020; Chen & Gnedin 2024) to make this connection. Having greater insight into the phase space of the distant debris of the trailing arm of the Sagittarius stream beyond the distances of NGC 2419 would yield useful information about the link between NGC 2419 and the Sagittarius stream. In general, prior work is limited to a subset of the stars in the distant trailing arm of the Sagittarius stream that only reach the apocentre of the trailing arm at most.

For this work, we used updated photometric catalogues and selection techniques (Thomas et al. 2018; Starkenburg et al. 2019) to select the best candidate BHB member stars of the outermost trailing arm apocentre and spur feature, and followed these stars up spectroscopically to determine their LOS velocities and constrain this feature kinematically for the first time. The most luminous standard candles observable with sufficient signal-to-noise ratio are outside of the limits of current spectroscopic surveys. Therefore, there is little published data on LOS velocities of bright stars in this outer realm of the Galaxy.

This work has two primary aims. Firstly, we investigate,
using spectroscopic data, the effectiveness of selecting BHB
stars using narrow- and/or broad-band photometry, especially at
the faint-end magnitudes, corresponding to BHB stars located
in the outermost region of the MW stellar halo at distances
above 50 kpc. To this end, we employed the ESO/VLT/FORS2
spectrograph with the 600B+22 grating covering the Balmer

<footnote n="1">¹ See Appendix in Belokurov et al. (2014) for this coordinate transformation.</footnote>
lines of the stars from which their discriminatory surface gravity
can be determined, separating BHB stars from contaminants<fnref n="2" />.
Secondly, we wanted to ascertain whether the LOS velocities
derived from the optical VLT/FORS2 spectra of the confirmed
BHB stars can provide further observational constraints on the
MW gravitational potential at large distances.

The structure of this study is as follows. In Sect. 2, we lay out
the characteristics of the sample of candidate BHB stars in the
Sagittarius stream spur feature, describe the methodology used
for the data reduction, spectral analysis, and stellar classification,
and present the data on the kinematics of the classified BHB
stars. We present and discuss the results in comparison to sim-
ulations and its implications for the Galactic potential at large
distances (Sect. 3.2) and to the MW globular cluster NGC 2419
(Sect. 3.3) in Sect. 3, and summarise in Sect. 4.

## 2. Data

### 2.1. Photometric selection of candidate blue horizontal branch stars

In order to cleanly select BHB stars, one needs to discrimi-
nate them from contamination sources with similar broad-band
colours, such as white dwarfs, quasars, and – most importantly
at the magnitudes of our study – blue straggler stars. Previous
research has established that there are various ways to deselect
spectroscopically and/or photometrically the mentioned contam-
inants (e.g. Xue et al. 2008, 2011; Vickers et al. 2012; Barbosa
et al. 2022; Yu et al. 2024; Amarante et al. 2024; Fukushima
et al. 2025; Byström et al. 2024).

Regarding quasars as contamination, analysis of the dered-
dened ($u-g$)$_{SDSS,0}$ versus ($g-r$)$_{SDSS,0}$ diagram of stellar-like
objects by Yanny et al. (2000) showed that A-type stars have
similar ($g-r$)$_{SDSS,0}$ colour index values to quasars, but are about
1 mag redder in ($u-g$)$_{SDSS,0}$. Using this or a similar colour space,
previous works have been able to deselect quasars as well as
white dwarfs (e.g. Xue et al. 2008, 2011; Deason et al. 2011,
2014; Ibata et al. 2017; Starkenburg et al. 2019; Fantin et al. 2019;
Barbosa et al. 2022).

Other works have looked at colour spaces using the near-
infrared instead of the near-ultraviolet in combination with opti-
cal photometry and report also success in filtering out quasars
and white dwarfs from (candidate) A-type stars (e.g. Vickers
et al. 2012; Deason et al. 2018; Fukushima et al. 2018, 2019; Yu
et al. 2024; Amarante et al. 2024; Fukushima et al. 2025). A key
study of Vickers et al. (2012) comparing SDSS ($i-z$)$_{SDSS,0}$ and
($g-r$)$_{SDSS,0}$ of stellar-like objects with spectroscopic data found
that quasars and white dwarfs are on average the reddest and
bluest, respectively, in ($i-z$)$_{SDSS,0}$ of all contaminants. To further
distinguish BHB stars from quasars, Vickers et al. (2012), in turn
based on the results of Lenz et al. (1998), also use ($g-r$)$_{SDSS,0}$,
since quasars tend to have higher values for this colour index
(especially at higher redshifts) due to more uniform emission in
these bands, in contrast to blue stars. Most studies investigating
the photometric selection of (candidate) BHB stars have been
based on such a combination of optical and near-infrared filters
using either data from the Subaru Telescope/Hyper Suprime-
Cam (e.g. Deason et al. 2018; Fukushima et al. 2018, 2019, 2025)
or the Dark Energy Survey (e.g. Yu et al. 2024).

To tackle the contamination from blue straggler stars at sim-
ilar magnitudes, three main ways are suggested in the literature.

Firstly, some studies have used typically low-resolution optical
spectra of A-type stars to use the measured widths of Balmer
lines as discriminator between BHB and more compact blue
straggler stars. Due to their higher surface gravity and thus
stronger pressure broadening, the Balmer lines are broader for
blue straggler stars than for BHB stars (see e.g. Clewley et al.
2002; Xue et al. 2008, 2011). A discrimination based on the
widths of Balmer lines can be expressed in these measured
widths directly, but also in resulting effective temperatures and
surface gravities estimated from (these same lines in the) spec-
tra (e.g. Barbosa et al. 2022; Byström et al. 2024). Vickers et al.
(2021) demonstrate that also a classifier trained on Large Sky
Area Multi-Object Fiber Spectroscopic Telescope (LAMOST,
Cui et al. 2012) spectra of BHB and non-BHB stars with suf-
ficient signal-to-noise ratio can yield a selection of (candidate)
BHB stars with a purity of 86 per cent.

Secondly, as suggested by Lenz et al. (1998) a near-ultraviolet
band as SDSS *u* that captures the Balmer break in hot stars can
also be used to differentiate between BHB and blue straggler
stars since the strength of the Balmer break in hot stars is depen-
dent on the surface gravity (see works by, e.g. Ruhland et al.
2011; Deason et al. 2011, 2014; Thomas et al. 2018; Starkenburg
et al. 2019).

Thirdly, it is now well established from a variety of studies
that also near-infrared filters around the Paschen lines in com-
bination with optical photometry also separates blue straggler
from BHB stars when the signal-to-noise ratio is high enough
(e.g. Vickers et al. 2012; Deason et al. 2018; Fukushima et al.
2018, 2019; Thomas et al. 2018; Yu et al. 2024; Amarante et al.
2024; Fukushima et al. 2025).

In summary, it has been shown that even photometric selec-
tions of candidate BHB stars can provide high purities between
80 and 90 per cent. The approach we follow in this work for pho-
tometric selection of candidate BHB stars in the spur feature of
the Sagittarius stream uses this technique as well.

While Thomas et al. (2018) focuses on the photometric selec-
tion of candidate BHB stars based on CFIS-$u_0$ and PS1 $g_0 r_0 i_0 z_0$,
Starkenburg et al. (2019) demonstrated the efficacy of employ-
ing a combination of broad-band SDSS $ugr_0$ and narrow-band
Pristine Survey CaHK$_0$ photometry (Starkenburg et al. 2017) to
identify candidate BHB stars in the Galactic halo and isolate
them from contaminating samples. The high efficiency of the
combination of the two blue *u* and CaHK bands is illustrated in
the lower middle and right panel of Fig. 1. Taken together, the *u*
and CaHK band demarcate the Balmer break and are therefore
highly effective to measure the higher surface gravity of blue
straggler compared to BHB stars. The lower middle and right
panel of Fig. 1 also indicate that the BHB and blue straggler stars
are closer to each other for the $u_{CFIS}$–Pristine CaHK colour index
than they are for $u_{SDSS}$–Pristine CaHK, as was already reported
by Titulaer (2021). The prime cause of this discrepancy is the
difference in wavelength range of the CFIS-*u* and SDSS-*u* filters
because CFIS-*u* is redder than SDSS-*u*.

### 2.2. Sample selection of candidate blue horizontal branch stars

Following the approach of Thomas et al. (2018), we initially
select 974 candidate BHB stars in the Galactic halo. These
candidate BHB stars were originally isolated from contaminat-
ing candidate blue straggler stars in a photometric sample that
combined PS1 *griz* data taken as part of the 3Π survey with
*u*-band data taken with CFHT/MegaCam. The latter observa-
tions were taken as part of the Canada-France Imaging Survey
(CFIS, Ibata et al. 2017), that has since been absorbed into the

<footnote n="2">² Based on observations collected at the European Southern Observa-
tory under ESO programmes 106.21L8.001 and 106.21L8.002 (PI: E.
Starkenburg).</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Overview of the target stars in multiple colour-colour spaces. The data were gathered from multiple sources: the Panoramic Survey Telescope and Rapid Response System 1 3Π survey Data Release 2 (Chambers et al. 2016; Magnier et al. 2020); the Canada-France Imaging Survey; the Sloan Digital Sky Survey Data Release 18 (Almeida et al. 2023); and the Pristine Survey (Starkenburg et al. 2017). The candidate BHB stars in our sample are colour-coded according to the ID number given in Table 1. The background stars come from two samples. The pink contours are blue straggler stars from the Xue et al. (2008) set while the brown contours represent BHB stars from the Xue et al. (2011) dataset (see the markers of the same colours in the upper left and upper right corners in the top and lower panels, respectively, for the average uncertainties of these stars). We use different markers to indicate the source of each star (Thomas et al. (2018), Starkenburg et al. (2019), or both catalogues). We decided to not show the stars with red colour indices ($u_{SDSS}$-Pristine CaHK)$_0$ in the lower middle panel that came from the catalogue Thomas et al. (2018) and just represent them with triangles pointing to the right.</figcaption></figure>

Ultraviolet Near Infrared Optical Northern Survey (UNIONS; Gwyn et al. 2025). In addition to CFIS, UNIONS consists of members of the Pan-STARRS team, and the Wide Imaging with Subaru HyperSuprime-Cam of the *Euclid* Sky (WISHES) team. CFHT/CFIS is obtaining deep *u* and *r* bands; Pan-STARRS is obtaining deep *i* and moderate-deep *z* band imaging, and Subaru is obtaining deep *z*-band imaging through WISHES and *g*-band imaging through the Waterloo-Hawaii IfA *g*-band Survey (WHIGS). These independent efforts are directed, in part, to securing optical imaging to complement the *Euclid* space mission (Euclid Collaboration: Mellier et al. 2025), although UNIONS is a separate collaboration aimed at maximising the science return of these large and deep surveys of the northern skies. We note that with the advent of these deeper bands using other facilities, the analysis of Thomas et al. (2018) could be redone using deeper data. However, here we use the existing BHB catalogue and demonstrate it is more than sufficient to meet our science goals.

The colour spaces, and the location of spectroscopically confirmed BHB and blue straggler stars (taken from the studies of Xue et al. 2008, 2011) are shown in the upper panels of Fig. 1. From the original 974 stars, 29 candidate BHB stars were

identified to be possibly associated to the spur feature due to their proximity in terms of declination (25–40 deg) and distance (~80–150 kpc) to this distant arm of the Sagittarius tidal stream. An additional 40 BHB candidates could be added in the same region from the sample of Starkenburg et al. (2019).

Figure 1 illustrates the power of the combined colour selection from photometry of DR2 of the PS1 3Π survey (Chambers et al. 2016), CFIS (Ibata et al. 2017), SDSS DR18 (Almeida et al. 2023), and the Pristine Survey (Starkenburg et al. 2017). For the estimation of interstellar extinction we used the integrated colour excess values along the LOS of each star from Schlegel et al. (1998) in combination with the tabulated values $A_X/E(B-V)_{SFD}$ from Schlafly & Finkbeiner (2011) assuming $R_V = 3.1$ and $A_{Pristine CaHK}/E(B-V)_{SFD} = 3.918$ ((Martin et al. 2024)), where we assume that extinction in CFIS-*u* is roughly the same as in Sloan-*u* following Ibata et al. (2017), Thomas et al. (2018, 2019, 2020), and Jensen et al. (2021).

Out of these 69 candidate BHB stars, 25 were followed up spectroscopically using FORS2. These were selected on the basis of their colour information. Of the 25 targets chosen, 17 had CaHK photometry from the A-star selection of Starkenburg et al. (2019), the other stars chosen were not present
<table><caption>**Table 1.** Basic properties of the candidate BHB stars in the sample.</caption><thead><tr><th>Designation</th><th>ID</th><th>RA<br>J2000</th><th>Dec<br>J2000</th><th>Gaia source ID</th><th>Observation time<br>Julian Date-UTC</th><th>Source</th></tr></thead><tbody><tr><td>144511105510814456</td><td>candBHB1</td><td>7h22m12.24s</td><td>30°25′40.8"</td><td>885923969299049856</td><td>2459145.820359</td><td>C</td></tr><tr><td>J072805.06+254349.2</td><td>candBHB2</td><td>7h28m05.04s</td><td>25°43′48"</td><td>...</td><td>2459193.744190</td><td>P</td></tr><tr><td>141301122397717196</td><td>candBHB3</td><td>7h28m57.6s</td><td>27°45′21.6"</td><td>873208357842184448</td><td>2459200.699201</td><td>P</td></tr><tr><td>147381140874742260</td><td>candBHB4</td><td>7h36m20.88s</td><td>32°49′04.8"</td><td>...</td><td>2459193.787593</td><td>C</td></tr><tr><td>J074315.57+331906.5</td><td>candBHB5</td><td>7h43m15.6s</td><td>33°19′04.8"</td><td>...</td><td>2459196.707963</td><td>P</td></tr><tr><td>J075310.18+271300.1</td><td>candBHB6</td><td>7h53m10.08s</td><td>27°13′01.2"</td><td>...</td><td>2459202.715949</td><td>P</td></tr><tr><td>148711185503191921</td><td>candBHB7</td><td>7h54m12s</td><td>33°55′33.6"</td><td>...</td><td>2459197.740104</td><td>C</td></tr><tr><td>J075724.50+330505.6</td><td>candBHB8</td><td>7h57m24.48s</td><td>33°05′06"</td><td>881458199742504832</td><td>2459225.660648</td><td>P</td></tr><tr><td>J075815.87+364537.2</td><td>candBHB9</td><td>7h58m15.84s</td><td>36°45′36"</td><td>906891484241931264</td><td>2459196.753993</td><td>PC</td></tr><tr><td>J075922.63+265611.1</td><td>candBHB10</td><td>7h59m22.56s</td><td>26°56′09.6"</td><td>874378478731783424</td><td>2459192.777477</td><td>P</td></tr><tr><td>J080051.99+361303.2</td><td>candBHB11</td><td>8h00m52.08s</td><td>36°13′04.8"</td><td>906624234198306816</td><td>2459200.766829</td><td>P</td></tr><tr><td>J080146.11+343120.7</td><td>candBHB12</td><td>8h01m46.08s</td><td>34°31′19.2"</td><td>905620380080032768</td><td>2459199.732419</td><td>P</td></tr><tr><td>J080152.47+253951.0</td><td>candBHB13</td><td>8h01m52.56s</td><td>25°39′50.4"</td><td>682026890911273088</td><td>2459176.742639</td><td>P</td></tr><tr><td>145911206860680542</td><td>candBHB14</td><td>8h02m44.64s</td><td>31°35′31.2"</td><td>...</td><td>2459198.756331</td><td>C</td></tr><tr><td>142131209582933555</td><td>candBHB15</td><td>8h03m49.92s</td><td>28°26′38.4"</td><td>876085642331737984</td><td>2459192.730613</td><td>C</td></tr><tr><td>J080408.41+355917.0</td><td>candBHB16</td><td>8h04m08.4s</td><td>35°59′16.8"</td><td>906925603461957248</td><td>2459202.755405</td><td>P</td></tr><tr><td>J080712.26+300447.1</td><td>candBHB17</td><td>8h07m12.24s</td><td>30°04′48"</td><td>876855197392497024</td><td>2459201.713171</td><td>P</td></tr><tr><td>J080757.78+350947.9</td><td>candBHB18</td><td>8h07m57.84s</td><td>35°09′46.8"</td><td>905889932227653888</td><td>2459224.701400</td><td>P</td></tr><tr><td>J080827.57+301532.4</td><td>candBHB19</td><td>8h08m27.6s</td><td>30°15′32.4"</td><td>876954290877415424</td><td>2459203.723565</td><td>PC</td></tr><tr><td>J080943.83+300352.5</td><td>candBHB20</td><td>8h09m43.92s</td><td>30°03′54"</td><td>876748269886574720</td><td>2459207.693669</td><td>PC</td></tr><tr><td>J081027.60+315155.4</td><td>candBHB21</td><td>8h10m27.6s</td><td>31°51′54"</td><td>901885613958625152</td><td>2459177.769769</td><td>P</td></tr><tr><td>J081507.72+294705.5</td><td>candBHB22</td><td>8h15m07.68s</td><td>29°47′06"</td><td>...</td><td>2459199.792697</td><td>P</td></tr><tr><td>142731243370082351</td><td>candBHB23</td><td>8h17m20.88s</td><td>28°56′34.8"</td><td>...</td><td>2459176.782975</td><td>C</td></tr><tr><td>J084934.60+375926.5</td><td>candBHB24</td><td>8h49m34.56s</td><td>37°59′27.6"</td><td>...</td><td>2459196.797465</td><td>P</td></tr><tr><td>152761351890263748</td><td>candBHB25</td><td>9h00m45.36s</td><td>37°18′10.8"</td><td>...</td><td>2459197.792465</td><td>P</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>2459195.808785</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>2459198.799965</td><td></td></tr><tr><td></td><td></td><td></td><td></td><td></td><td>2459225.707431</td><td>C</td></tr></tbody></table>

**Notes.** We provide the data on equatorial coordinates in the International Celestial Reference System from PS1 DR2 (Chambers et al. 2016; Magnier et al. 2020). We use the SDSS and PS1 designation in the first column if the star came from the Pristine Survey ('P' in last column) (Starkenburg et al. 2019), the CFIS catalogue ('C' in last column) (Thomas et al. 2018), or from both catalogues ('PC').

in this photometric catalogue at the time. While 24 stars have CFIS-u magnitudes, 23 stars have SDSS-u magnitudes. For the remainder of the paper, we refer to them with their internal IDs for this study (second column in Table 1).

Figure 1 shows the location of the 25 chosen targets in the various colour-colour spaces that were used for their selection on top of the spectroscopically identified samples. Because of the pilot-programme nature of this observational campaign, it was not enforced that the stars were the best BHB-candidates in all possible colour spaces. The programme was deliberately leaving the option open that one of the colour spaces would deliver an inferior classification. Additionally, we note that updates in the PS1, Pristine survey, and CFIS-u more recent (internal) data releases are affecting the colours with respect to the original selection. For Pristine the main updates have been an increase of footprint, along with a re-calibration of the narrow-band CaHK photometry and updates to the data reduction pipeline (Martin et al. 2024) and for CFIS the changes between the original catalogue used by Thomas et al. (2018) and the latest data reduction are due to a change of the photometric reduction pipeline and a better estimation of the zero point. The changes in colour for the targets and background of (candidate) A-type stars in Xue et al. (2008) and Xue et al. (2011) due to these updates are typically small, as can be appreciated from the fact that in Fig. 1 still most targets lie in the regions covered by BHB stars.

These updates have added some target stars from the original catalogue by Thomas et al. (2018) to the lower middle panel of Fig. 1 (e. g., candBHB3, candBHB7, and candBHB25; represented by circles) and some programme stars from the Pristine Survey catalogue (Starkenburg et al. 2019) to the upper panels of Fig. 1 (indicated by octagons) although these colours/colour-colour-diagrams were originally not used for the selection of these stars. Most strikingly, we see that candBHB25 is very discrepant in the colour-colour diagrams of the lower and upper left panels compared to where (candidate) BHB stars are to be expected in these projected colour spaces. By contrast, we can also see that the same star is significantly closer to the sequence of (candidate) BHB stars in the upper right panel of Fig. 1. Moreover, candBHB24 and candBHB5 are outliers in the top left panel, but fall nicely on the sequence of candidate BHB stars in several other panels.

We note that the majority of candidate BHB stars in the sample have very faint PS1 g magnitudes for their Data Release 2, going well beyond 20 (up to 21.6 AB mag). A consequence of this faint selection is the limited overlap of the target stars in Gaia DR 3, due to incompleteness beyond Gaia G ≈ 21 mag (Gaia Collaboration 2023) (see also Sect. 3.3). The detected, candidate BHB stars in Gaia DR3 (Gaia Collaboration 2023) have all parallaxes consistent with being unconstrained/unresolved. This is further evidence that
<table><caption>**Table 2.** Key characteristics of the velocity standard stars.</caption><thead><tr><th>Designation</th><th>RA<br>J2016</th><th>Dec</th><th>Spectral type</th><th>Barycentric LOS velocity<br>[km s<sup>-1</sup>]</th></tr></thead><tbody><tr><td>BD +24 1843</td><td>8h04m42.70936675s</td><td>24°19′49.47351216″</td><td>G</td><td>25.168 ±0.004</td></tr><tr><td>BD +31 1781</td><td>8h18m10.42693323s</td><td>30°35′49.7809604″</td><td>K</td><td>13.457 ±0.019</td></tr><tr><td>BD +34 1955</td><td>9h12m37.57008557s</td><td>33°36′01.09932489″</td><td>K</td><td>1.539 ±0.004</td></tr><tr><td>WDS J07277+2420A</td><td>7h27m39.94289657s</td><td>24°20′09.97879538″</td><td>K</td><td>-18.559 ±0.208</td></tr><tr><td>BD +26 1647</td><td>7h46m58.48254722s</td><td>26°01′28.79033441″</td><td>G</td><td>13.109 ±0.008</td></tr><tr><td>BD +31 1684</td><td>7h53m33.99666092s</td><td>30°35′48.89008376″</td><td>G</td><td>-234.195 ±0.013</td></tr><tr><td>BD +30 1501</td><td>7h26m19.81197854s</td><td>29°58′07.98205929″</td><td>K</td><td>-36.684 ±0.044</td></tr><tr><td>TYC 2461-988-1</td><td>7h32m44.21719125s</td><td>33°50′06.07594227″</td><td>F</td><td>23.187 ±0.020</td></tr><tr><td>BD +36 1823</td><td>8h27m18.62972372s</td><td>35°49′01.65250699″</td><td>F</td><td>9.282 ±0.043</td></tr><tr><td>BD +35 1801</td><td>8h19m27.87619964s</td><td>35°01′21.76871662″</td><td>F</td><td>-23.908 ±0.013</td></tr></tbody></table>

Notes. The data were sourced from Gaia Collaboration (2023) except the barycentric LOS velocities that are from Soubiran et al. (2018). We used the same approach for the designations as noted in Table 1. The equatorial coordinates are also in the International Celestial Reference System. Spectral types are from the Extended Stellar Parametrizer for Hot Stars module as part of Gaia DR3 Astrophysical parameters inference system (Creevey et al. 2023).

these stars are at (large) distances where Gaia is not able to constrain/resolve these parallaxes.

## 2.3. VLT/FORS2 spectroscopic set-up

We collected VLT/FORS2 data in the wavelength range 3300–6210 Å using the low spectral resolution ($R = 780$ at 4627 Å with a 1 arcsec) 600B+22 grism for the candidate BHB stars listed in Table 1. The slit width was 1 arcsec, exposure times are close to 2300 s, and observational epochs are listed in Table 1 (Period 106, PI: E. Starkenburg; ESO programmes 106.21L8.001 and 106.21L8.002). We note that Deason et al. (2012) state that uncertainties in the LOS velocities on the level of $6 \text{ km s}^{-1}$ can be reached with the FORS2 instrumental set-up we use for this work.

A key ingredient for reliable and accurate LOS velocity determinations is the quality of the wavelength calibration. In the FORS2 spectrograph, the wavelength calibration is based on the lines in spectra of arc-lamps taken at daytime and at zenith. However, for spectra of stars with significant zenith distance (such as at high air masses as in our sample), there is a potential for bias in the wavelength calibrations due to instrument flexure under gravity (Anderson 2020). Strong levels of instrument flexure can give rise to second-order effects in the wavelength solution.

Many researchers have utilised night sky emission lines and/or telluric lines to measure these second-order effects in the dispersion solution that are gathered together with the science spectra in the same set-up and under the same influence of instrument flexure (e.g. Southworth et al. 2006; Yanny et al. 2009; Deason et al. 2012; Caffau et al. 2020). However, in the wavelength range of our study, only a few suitable lines are available. Deason et al. (2012) investigated the differential impact of these second-order effects using strong night sky emission lines in FORS2 spectra of 48 candidate BHB stars in a set-up similar to that used in this programme. A more recent study by Caffau et al. (2020) also involved a similar FORS2 set-up to this programme's, and reported more variation of the LOS velocities of their stars with reference values, after a shift to the dispersion axis of their FORS2 spectra based on the measured positions of eight unblended night sky emission lines. Following this result, we refrain from this approach.

Instead, to establish whether instrument flexure had a significant impact on the wavelength solution derived from arc-lamps

taken at daytime and at zenith, we additionally gathered data with the same spectral set-up (except a smaller slit of 0.4 arcsec) for ten standard stars. Each observing block of our observing programme would typically have a BHB candidate target, immediately followed or preceded by the observation of a velocity standard star as close in the sky as possible. The list of velocity standards, and their LOS velocity measurements from literature, are presented in Table 2 and taken from Soubiran et al. (2013, 2018). This approach provides us with a second avenue to test the accuracy and precision of the measured LOS velocities.

## 2.4. Spectroscopic reduction

Standard two-dimensional, CCD data reduction and calibration were performed using the FORS2 pipeline 5.5.7 (2021) within the ESOREFLEX 2.11.3 (2021) environment (Freudling et al. 2013). Details are given in Appendix A<fnref n="3" />.

The results of the reduction and calibration are set out in Fig. 2 for three example spectra of candidate BHB stars in our sample. While most of the targets were observed with one single exposure, our faintest programme targets have multiple exposures. The spectra at 4500 Å have signal-to-noise ratios in the range 7–25. This figure depicts the general pattern of a prominent blue continuum and Balmer lines of most of the shown spectra of the candidate BHB stars, as we would expect for these types of stars (see some example lines of the Balmer series in the right panels).

The strong and narrow peaks at various wavelengths in the top two shown spectra in Fig. 2 are caused by pixels in the CCD data in the area of the spectra affected by cosmic rays. While it would have been an option to include a removal of the effects of cosmic rays in the FORS2 pipeline in a fore- and background subtraction based on a global sky spectrum, a global sky subtraction does not take into account variability of the spectral resolution across the slit axis. In addition, given that the cosmic rays present in the spectra of some candidate BHB stars do not affect all the Balmer lines that will be used to identify true

<footnote n="3">³ Data of the first observation of candBHB1 was sourced from the ESO archive later as part of the open stream release FOR S2-SPEC. Based on data obtained from the ESO Science Archive Facility with DOI: https://doi.eso.org/10.18727/archive/77</footnote>
<figure><img src="image_4.png" /></figure>

**Fig. 2.** Example spectra of three candidate BHB stars listed in Table 1 are shown in the left panels where the wavelength scale is relative to the spectrograph's reference frame at the time of each observation. The background shaded areas represent the uncertainties of the spectra. The solid red lines indicate the best-fitting models. The labels at the top annotate some of the strong Balmer lines visible in the spectra that can also be seen in the right panels for each example spectrum. We show the spectrum of the star with the lowest, median of, and highest median signal-to-noise ratio of all spectra across the dispersion axis considered for the fitting in the top, middle, and bottom panel, respectively.

BHB stars and misclassified BHB stars, we choose to keep them in the spectra.

The calibrated fluxes rise and fall steeply around 5570–5580 Å for the shown spectra in Fig. 2. This is probably due to the non-optimal performance of the sky subtraction in that wavelength region where the strong [O I] night sky emission line is located. We note here again that this does not affect the lines of the Balmer series. That is why we left this feature in the spectra as it is.

In Fig. 2, there is a clear trend of decreasing signal-to-noise noise blueward of around 3800 Å. This could be attributed to the following factors. Firstly, it is known that the blue flat field lamp is unstable below 3800 Å leading to variations of the spectral energy distribution of 20 per cent and systematic distortions in the calibrated flux given that flat fields are obtained separately for science and flux standard star observations (Anderson 2020). Secondly, extrapolation of the wavelength solution for the blue end due to a lack of spectral lines in the arc-lamps at the blue edge results in significant residuals in this part of the spectrum. This does not affect the central regions of the spectra where most of the Balmer lines are located. It can therefore be assumed that the noise-dominated, blue part of the spectrum has no impact on our analysis. In the remainder of this work, we considered only 3900 Å to 5990 Å to be the reliable range of the spectrum.

## 2.5. Spectral fitting

We used RVSPECFIT (Koposov et al. 2011; Koposov 2019)<fnref n="4" />, to analyse the spectra. A major advantage of spectral fitting with

RVSPECFIT is that it provides both estimates of the LOS velocity and of effective temperature and surface gravity from each spectrum. Whereas the first parameter provides us the data to constrain the velocity structure of the spur feature of the Sagittarius stream, the latter makes it possible to classify our sample of candidate BHB stars into bona fide BHB, blue straggler stars, or other contaminants in a way similar to that done by Barbosa et al. (2022) and Byström et al. (2024) (the latter also uses results from RVSPECFIT).

To estimate the LOS velocity, effective temperature, and surface gravity of a star from a spectrum, RVSPECFIT determines the best-fitting synthetic spectrum from a set of interpolated synthetic stellar spectra from the PHOENIX library (Husser et al. 2013), version 2. The shift results in a measured apparent LOS velocity of the star. Any remaining imprecision in the data is taken into account by a multiplicative, normalising polynomial. In greater detail, RVSPECFIT maximises the multivariate Gaussian likelihood that describes a stellar spectrum. RVSPECFIT utilises the effective temperature $T_{\text{eff}}$, $\log_{10}(g)$, [Fe/H], and [α/Fe] as parameters of the PHOENIX library of templates and LOS velocity<fnref n="5" />. It arrives at these results after marginalising over the coefficients of the normalising polynomial, as presented in

<footnote n="5">⁵ While RVSPECFIT includes the capability to estimate the stellar rotation parameter $v \sin i$, in section 6 about the known limitations of the Dark Energy Spectroscopic Instrument Early DR based also on results from RVSPECFIT Koposov et al. (2024) recommend not to use estimates $v \sin i$ from the current version of RVSPECFIT due to several issues as mentioned in Koposov et al. (2024). Because of this and also the irrelevance of this parameter for the scope of this work, we do not consider the output $v \sin i$ by RVSPECFIT.</footnote>

<footnote n="4">⁴ Version 0.4.0.240116+dev at github.com/segasai/rvspecfit</footnote>
the initial $\chi^2$ of the spectroscopic data, given a template of the prepared and interpolated PHOENIX grid. Ju et al. (2024) studied the effects of varying $T_{\text{eff}}$, [Fe/H], and signal-to-noise ratio S/N on typical spectra of BHB stars and note that the Calcium II K line is the only metallicity-sensitive line that has significant strength in the optical range. Additionally, they show that even this strong metal line gets very weak for typical BHB stellar parameters, as is evidenced by the fact that for a BHB star with $[\text{Fe/H}] = -1.7$ and $T_{\text{eff}} = 8000$ K the line is not visible any more for S/N $\le 20$ and low spectral resolution (see their Figure 14). For slightly hotter stars, with $T_{\text{eff}} = 9000$ K, this line is not detectable at similar conditions already at $[\text{Fe/H}] = -1$ (see Ju et al. 2024, Figure 15). These examples illustrate that it is not possible to robustly measure $[\text{Fe/H}]$ and $[\alpha/\text{Fe}]$ for spectra of the resolution and signal-to-noise in our sample. In order to quantify the robustness of the resulting $T_{\text{eff}}$ and $\log_{10}(g)$ from RVSPECFIT to different values of $[\text{Fe/H}]$ and $[\alpha/\text{Fe}]$ we run the spectral fitting for the assumed, fixed values of $[\text{Fe/H}] = -3.8$ (one step size away from the edge grid), $[\text{Fe/H}] = -1.8$ (close to the middle of the grid), and $[\text{Fe/H}] = 0.0$ and $[\alpha/\text{Fe}] = 0.0$ and found no significant differences between the resulting effective temperature and surface gravity from the fits for each assumed, fixed value of $[\text{Fe/H}]$ and $[\alpha/\text{Fe}]$.

To begin the process of spectral fitting with RVSPECFIT, we had to prepare the set of synthetic stellar spectra from the PHOENIX library, in the three instrumental set-ups that we use and interpolate between the grid parameters. The three instrumental set-ups that we use are connected to the spectra of the candidate BHB stars in our sample taken with the slit with a width of 1.0 arcsec, the spectra of the velocity standard stars taken with the 0.4 arcsec slit, and Sloan Digital Sky Survey DR17 spectra (Abdurro'uf et al. 2022) of a set of reference BHB and blue straggler stars that we use to classify our target stars in the effective temperature versus surface gravity space (Kiel diagram).

The process of parameter estimation was carried out given the observed spectra and priors for each parameters. By default, RVSPECFIT assumes uniform priors of $T_{\text{eff}}$, $\log_{10}(g)$, [Fe/H], and $[\alpha/\text{Fe}]$ over the entire PHOENIX grid although Gaussian priors can be custom defined. As detailed below, we use this option for parameters we have prior information for. We additionally adopt a broad uniform prior of the LOS velocity between $-500 \text{ km s}^{-1}$ and $500 \text{ km s}^{-1}$ throughout all performed fitting with RVSPECFIT to make no strong prior assumptions for the most important parameter of interest in this study.

For the spectra of the velocity standard stars, we determined LOS velocities with RVSPECFIT using known stellar parameters ($T_{\text{eff}}$, $\log_{10}(g)$, [Fe/H], $[\alpha/\text{Fe}]$) as priors. More details on this will be given in the appendix (Appendix A).

As shown in Table 3, we do not set priors on $\log_{10}(g)$ for our programme stars as we measure this value to make the BHB-blue straggler distinction. We have information on $T_{\text{eff}}$ from their colours (see Appendix B) and adopt a (conservative) prior based on this information.

Finally, we set the priors included in the posterior calculation when fitting the SDSS spectra of BHB and (candidate) blue straggler stars in the Barbosa et al. (2022) and Xue et al. (2008) catalogues, respectively, with RVSPECFIT. The normally distributed priors used here were synthesised using the same method detailed for the spectra of the candidate BHB stars, with effective temperature values sourced from two sources, depending on whether it is a BHB or (candidate) blue straggler star, from the Barbosa et al. (2022) or Xue et al. (2008) sets, respectively. These are the results of the Sloan Extension

of Galactic Understanding and Exploration (SEGUE) Stellar Parameter Pipeline, as presented by Lee et al. (2008a,b) and based on SEGUE data (Yanny et al. 2009) for the BHB stars. In the Xue et al. (2008) sample, we adopt the effective temperature of the model star with the ELODIE template Prugniel & Soubiran (2001) that best describes the SDSS spectrum of each (candidate) blue straggler star. The latter values are stored as metadata in the SDSS spectral data. In contrast to the former source, the metadata only have the effective temperatures without uncertainties. Consequently, we assumed a uniform, conservative standard deviation of 1000 K for all normally distributed priors of effective temperature when fitting SDSS spectra of the (candidate) blue straggler stars with RVSPECFIT.

Another significant aspect of the spectral fitting with RVSPECFIT is the set of radial basis functions used to describe the normalising polynomial. A set of ten radial basis functions was used, following the procedure described in Cooper et al. (2023) used by the MW survey pipeline of the Dark Energy Spectroscopic Instrument (DESI-MWS).

RVSPECFIT for all our targets outputs a maximum a posteriori estimation of the LOS velocity while also computing higher moments of the projection of the posterior across the velocity axis (including standard deviation). The estimates of the uncertainty of the maximum a posteriori values of $T_{\text{eff}}$ and $\log_{10}(g)$ show standard deviations of a Gaussian that was used to approximate the four-dimensional posterior of ($T_{\text{eff}}$, $\log_{10}(g)$, [Fe/H], $[\alpha/\text{Fe}]$) around the mode. This approximation is equivalent to evaluating the Hessian matrix of the posterior at the mode to obtain the covariance matrix of the Gaussian (Bailer-Jones 2017) as done in RVSPECFIT.

## 2.6. Line-of-sight velocities and assessment of accuracy of their measurement

As our data reduction does not include any heliocentric correction and conversion to the Local Standard of Rest, they are in the observer's reference frame. Heliocentric velocity corrections and conversion to the Local Standard of Rest were calculated using NOAO.RV.RVCORRECT in the Image Reduction and Analysis Facility (IRAF, Tody 1986, 1993; National Optical Astronomy Observatories 1999). All these calculations were carried out using PYRAF, that is, a command language for IRAF based on PYTHON. The corrections were computed using the coordinates of the stars in Tables 1 and 2, mid-observation times calculated from the data in Table 1, and total magnitude of the Cartesian, Galactic velocity (11.1, 12.24, 7.25) km s$^{-1}$ of the Sun with respect to the Local Standard of Rest as found by Schönrich et al. (2010).

The main source of error in our heliocentric LOS velocity estimates is a potential for bias resulting from the telescope's exposure to flexure under gravity during observation of high air mass target stars with the FORS2 instrument. This bias is not accounted for in the wavelength solution of the spectra, which is derived based on reference spectra of arc lamps obtained with the telescope pointing towards zenith. In order to assess whether and how the maximum a posteriori estimates of the heliocentric LOS velocities (hereafter LOS velocities) are affected, we measured LOS velocities from the fitting of the spectra of the velocity standard stars listed in Table 2 with RVSPECFIT and compared them to the fiducial values in Table 2 from Soubiran et al. (2018). As shown in Fig. 3, the observed LOS velocities of the velocity standard stars here are offset from the data given in the catalogue of Soubiran et al. (2018) on average by $10-25 \text{ km s}^{-1}$.
<table><caption>**Table 3.** Spectroscopic (effective temperature, log₁₀ g, LOS velocity) obtained with RVSPECFIT and photometric parameters (heliocentric distance) of the candidate BHB stars.</caption><thead><tr><th>ID</th><th>Observation time<br>Julian Date-UTC</th><th>T<sub>eff</sub><br>Prior<br>[K]</th><th>Result<br>[K]</th><th>log<sub>10</sub>(g)</th><th>V<sub>h</sub><br>[km s<sup>-1</sup>]</th><th>Class</th><th>D<sub>☉</sub><br>[kpc]</th></tr></thead><tbody><tr><td>candBHB1</td><td>2459145.820359<br>2459193.744190</td><td>8423 ±562</td><td>8245 ±253<br>8611 ±213</td><td>3.080 ±0.132<br>3.005 ±0.104</td><td>10 ±14<br>35 ±12</td><td>BHB</td><td>106 ±5</td></tr><tr><td>candBHB2</td><td>2459200.699201</td><td>8376 ±842</td><td>8611 ±219</td><td>3.229 ±0.108</td><td>35 ±11</td><td>BHB</td><td>123 ±7</td></tr><tr><td>candBHB3</td><td>2459177.726597</td><td>8140 ±594</td><td>7691 ±121</td><td>3.649 ±0.405</td><td>78 ±20</td><td>Straggler</td><td>nan ±nan</td></tr><tr><td>candBHB4</td><td>2459193.787593</td><td>8361 ±1362</td><td>8361 ±238</td><td>2.947 ±0.104</td><td>-31 ±14</td><td>BHB</td><td>123 ±8</td></tr><tr><td>candBHB5</td><td>2459196.707963</td><td>8317 ±1447</td><td>7628 ±85</td><td>2.766 ±0.133</td><td>50 ±12</td><td>BHB</td><td>135 ±7</td></tr><tr><td>candBHB6</td><td>2459202.715949</td><td>8354 ±994</td><td>7957 ±110</td><td>3.284 ±0.403</td><td>54 ±17</td><td>BHB?</td><td>141 ±7</td></tr><tr><td>candBHB7</td><td>2459197.740104</td><td>8144 ±854</td><td>8341 ±221</td><td>3.091 ±0.109</td><td>39 ±12</td><td>BHB</td><td>116 ±6</td></tr><tr><td>candBHB8</td><td>2459225.660648</td><td>8002 ±429</td><td>7737 ±59</td><td>3.007 ±0.195</td><td>27 ±9</td><td>BHB</td><td>94 ±4</td></tr><tr><td>candBHB9</td><td>2459196.753993</td><td>8624 ±773</td><td>9068 ±278</td><td>3.340 ±0.129</td><td>122 ±9</td><td>BHB</td><td>90 ±4</td></tr><tr><td>candBHB10</td><td>2459192.777477</td><td>8192 ±463</td><td>7828 ±82</td><td>2.852 ±0.083</td><td>136 ±8</td><td>BHB</td><td>95 ±5</td></tr><tr><td>candBHB11</td><td>2459200.766829</td><td>8009 ±889</td><td>8975 ±240</td><td>3.255 ±0.128</td><td>31 ±9</td><td>BHB</td><td>92 ±5</td></tr><tr><td>candBHB12</td><td>2459199.732419</td><td>7729 ±1129</td><td>7843 ±92</td><td>2.973 ±0.112</td><td>37 ±12</td><td>BHB</td><td>128 ±6</td></tr><tr><td>candBHB13</td><td>2459176.742639</td><td>8119 ±401</td><td>8121 ±121</td><td>2.990 ±0.077</td><td>9 ±9</td><td>BHB</td><td>88 ±4</td></tr><tr><td>candBHB14</td><td>2459198.756331</td><td>8006 ±975</td><td>8553 ±206</td><td>3.267 ±0.136</td><td>71 ±12</td><td>BHB</td><td>148 ±8</td></tr><tr><td>candBHB15</td><td>2459192.730613</td><td>8030 ±655</td><td>8467 ±52</td><td>4.125 ±0.109</td><td>-48 ±8</td><td>Straggler</td><td>nan ±nan</td></tr><tr><td>candBHB16</td><td>2459202.755405</td><td>8502 ±421</td><td>9007 ±256</td><td>3.216 ±0.134</td><td>43 ±12</td><td>BHB</td><td>96 ±5</td></tr><tr><td>candBHB17</td><td>2459201.713171</td><td>8817 ±440</td><td>8616 ±180</td><td>3.221 ±0.042</td><td>45 ±13</td><td>BHB</td><td>97 ±5</td></tr><tr><td>candBHB18</td><td>2459224.701400</td><td>8420 ±856</td><td>8158 ±98</td><td>3.200 ±0.089</td><td>87 ±8</td><td>BHB</td><td>97 ±5</td></tr><tr><td>candBHB19</td><td>2459203.723565</td><td>8072 ±676</td><td>8104 ±99</td><td>3.226 ±0.108</td><td>16 ±10</td><td>BHB</td><td>89 ±4</td></tr><tr><td>candBHB20</td><td>2459207.693669</td><td>7878 ±281</td><td>7760 ±83</td><td>5.089 ±0.216</td><td>57 ±31</td><td>Straggler</td><td>nan ±nan</td></tr><tr><td>candBHB21</td><td>2459177.769769</td><td>8140 ±704</td><td>8350 ±159</td><td>3.063 ±0.091</td><td>68 ±9</td><td>BHB</td><td>90 ±4</td></tr><tr><td>candBHB22</td><td>2459199.792697</td><td>7821 ±627</td><td>7513 ±59</td><td>2.855 ±0.270</td><td>54 ±12</td><td>BHB</td><td>140 ±7</td></tr><tr><td>candBHB23</td><td>2459176.782975</td><td>8567 ±1097</td><td>8752 ±354</td><td>3.027 ±0.161</td><td>55 ±12</td><td>BHB</td><td>129 ±8</td></tr><tr><td>candBHB24</td><td>2459197.792465<br>2459196.797465<br>2459195.808785</td><td>9756 ±1409</td><td>10415 ±852<br>8583 ±262<br>8935 ±358</td><td>3.873 ±0.367<br>3.380 ±0.105<br>3.329 ±0.211</td><td>76 ±15<br>70 ±14<br>73 ±15</td><td>BHB</td><td>150 ±8</td></tr><tr><td>candBHB25</td><td>2459225.707431<br>2459198.799965</td><td>8735 ±1372</td><td>10787 ±nan<br>11118 ±2</td><td>6.456 ±nan<br>6.465 ±0.002</td><td>100 ±28<br>124 ±29</td><td>WD</td><td>nan ±nan<br>nan ±nan</td></tr></tbody></table>

**Notes.** A small subset of the returned parameters was chosen because of the expected difficulty in obtaining reliable estimates of [Fe/H], [α/Fe], and v sin i from low-resolution optical spectra of hot stars with strong Balmer lines that blend some of the metallicity-sensitive lines (see also noted limitations of RVSPECFIT in Koposov et al. 2024). We report the heliocentric LOS velocities without any corrections based on the velocity standard stars.

Of interest here is the significant increase in the discrepancy between the estimates of the LOS velocities of the velocity standard stars and the fiducial values with decreasing air mass with an estimated Spearman's rank correlation coefficient $ρ = 0.44 \pm 0.03$ using the SCIPY.STATS.SPEARMANR function (Virtanen et al. 2020) and the Pearson's correlation coefficient of $0.34 \pm 0.03$ and slope of the linear regression of $28 \pm 2$ km s⁻¹ that is both output by the SCIPY.STATS.LINREGRESS function (Virtanen et al. 2020). The uncertainties given for all these three values are their Monte Carlo standard deviations derived from 5000 reruns of each correlation test based on the difference of the line-of-sight velocities and their uncertainties. This result is somewhat counterintuitive, as the expectation would be that with rising air mass stronger possible second-order residuals would affect the wavelength calibration of the spectra more as a result of variations caused by telescope flexure under gravity.

A possible explanation for the behaviour in Fig. 3 described above is that it is a linear trend superimposed on a systematic offset of the order of $10–25$ km s⁻¹ (mean and standard deviation of $-15.6$ km s⁻¹ and $11.7$ km s⁻¹, respectively).

A $10–25$ km s⁻¹ offset in this study corroborates earlier findings by Caffau et al. (2020). Caffau et al. (2020) found that as measured LOS velocities for four stars from data in a similar instrumental set-up to that used in our work are compared to reference (fiducial) values they are on average offset by $17$ km s⁻¹ with a standard deviation $= 15$ km s⁻¹. However, other studies suggested that such an effect could be also lower (around $6$ km s⁻¹, Deason et al. 2012).

Quantification of a trend with air mass must be approached with some caution because the strength of the variations of the dispersion axis of the spectra against some reference depends also on how the telescope was positioned before the observations of the velocity standard stars, that is, how much zenith distance it moved when targeting the stars. It is beyond the scope of this paper to examine in detail the effect of telescope motion across the sky (between observations) on our velocity measurements.

We investigated the relation between our derived LOS velocity and the reference value with observing time and S/N but found no clear trend (see Appendix C). Because a velocity standard is observed for each programme star, we use for our
<figure><img src="image_5.png" /></figure>

**Fig. 3.** Results obtained from the analysis of the FORS2 spectra of the velocity standard stars collected from the data of the observations in connection with the entries in Table 2 compared to the reference barycentric LOS velocity values by Soubiran et al. (2018) (also given in Table 2). We plot the estimates of the difference of heliocentric or barycentric LOS velocities against the air mass of the star at the time of observation.

results both the uncorrected LOS velocity as well as corrected LOS velocity by the offset of the velocity standard closest in time to each BHB star to investigate whether the LOS velocity correction significantly affects the results.

## 2.7. Classification of candidate blue horizontal branch stars in sample

We classify the candidate BHB stars in our sample into bona fida BHB stars and other contaminants by comparing the positions of the target stars in a Kiel diagram. To establish the relevant parameter space we reference fitting results of SDSS DR17 spectra of a set of BHB and (candidate) blue straggler stars from the Barbosa et al. (2022) and Xue et al. (2008) catalogues through the same method (RVSPECFIT). A major advantage of using effective temperature and surface gravity inferred from spectral fitting in contrast to utilising the Balmer line widths (see e.g. Clewley et al. 2002; Xue et al. 2008, 2011) is that it encapsulates all the information from the spectrum instead of only individual lines. Barbosa et al. (2022) and Byström et al. (2024) identified BHB stars in SEGUE (Yanny et al. 2009) and DESI (DESI Collaboration 2016), respectively, in a similar way.

Figure 4 shows reference results versus our sample of target stars in a Kiel diagram. We highlight the low number of contaminants in our sample. Of the study population, 21 target stars have (exposure-averaged, maximum a posteriori) estimates of $T_{\text{eff}}$ and $\log_{10}(g)$ consistent with the trend in Fig. 4 of spectroscopically classified BHB stars. The results indicate that candBHB6 could be also in the area of the Kiel diagram where we expect blue straggler stars given its relatively large uncertainties.

Table 3 shows the proportion of different categories of stellar classification of the target stars according to their positions in a Kiel diagram. The bona fide BHB stars are indicated by ‘BHB’ in the ‘Class’ column in Table 3. The contaminants in our sample classified as blue straggler stars (candBHB15, candBHB20, and

<figure><img src="image_6.png" /></figure>

**Fig. 4.** Trends identified in one projection of the stellar parameter space for the group of target stars in the spur feature of the Sagittarius stream with FORS2 spectra. A large majority of the estimations of effective temperature and surface gravity for the target stars from RVSPECFIT (Koposov et al. 2011; Koposov 2019) (open star markers show estimates from fitting spectra obtained from each exposure for candBHB1 and candBHB24 for which several exposures exist) are in the range expected from bona fide BHB stars (see circles). A few outliers are marked with the numbers and more closely match the cross symbols that represent a reference dataset of blue straggler stars in the Xue et al. 2008 sample. Another outlier, candBHB25, is a likely white dwarf candidate highlighted just with an arrow because the best-fitting $T_{\text{eff}}$ and $\log_{10} g$ for both spectra are close to the upper edge of the PHOENIX grid, and are therefore not reliable. The grey marker in the upper left corner provides an estimate of the average uncertainties of the reference dataset.

candBHB3 although with less certainty due to proximity to the sequence of BHB stars in Fig. 4) are marked with ‘Blue straggler’ in the last column in Table 3. One candidate is probably a white dwarf (candBHB25). It has unreliable RVSPECFIT output from the spectra gathered from the CCD data of two exposures at different dates because they are close to the upper limit of the PHOENIX grid (labeled with ‘WD’ and represented with an arrow in Fig. 4). Roughly a third of those candidate BHB stars from the Thomas et al. (2018) catalogue (36 per cent or four out of the eleven target stars listed in Table 1 with ‘C’ or ‘PC’) are contaminants. This is slightly larger than their estimated contamination fraction of 24 per cent, but this can be attributed to the low-number statistics of our sample. For the subset of 17 candidate BHB stars in our sample selected from the Starkenburg et al. (2019) catalogue, only 6 per cent are flagged as contamination.

The findings indicate that the purity of the sample in question (80 per cent) corresponds to the anticipated range for the photometric selection of candidate BHB stars, as outlined in Starkenburg et al. (2019). This is evident in the middle panel of Figure 5 in Starkenburg et al. (2019), which assumes a one-to-one ratio of BHB to blue straggler stars in the outer halo.

## 2.8. Blue horizontal branch stars in Sagittarius spur region

In the following, we present the kinematics of the stars classified as BHB stars. Figure 5 presents a map of the Sagittarius
<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. Sagittarius stream according to a sample of candidate members in sample of (candidate) variable RR Lyrae stars by Hernitschek et al. (2017) in one projection across longitude of the Sagittarius stream coordinate system as defined by Vasiliev et al. (2021) (implemented in GALA Price-Whelan 2017), and heliocentric distance. For each (candidate) RR Lyrae, we colour-code the marker by the latitude in the same Sagittarius stream coordinate system. The top half of the mapped Sagittarius stream at between 135° and 180° longitude (see lower panel for enlarged plot of this part), includes the BHB stars in our sample plotted with distinct colours based on their ID (in the label of the colour scheme we only show the number at the end of the ID). The grey-shaded regions (dashed lines) represent the regions where the sample of (candidate) variable RR Lyrae stars by Hernitschek et al. (2017) are incomplete.</figcaption></figure>

stream according to a sample of candidate members in sample of variable RR Lyrae stars by Hernitschek et al. (2017) in one projection across longitude of Sagittarius stream coordinate system as defined by Vasiliev et al. (2021) and heliocentric distance. We also show the BHB stars in the Sagittarius spur region in Fig. 5. To estimate the heliocentric distances of the BHB stars, numerous studies of the Galactic stellar halo traced by BHB stars have utilised the fourth-order polynomial absolute magnitude

$M_{g_{SDSS,0},BHB}$ relation derived from the SDSS DR6 photometry of BHB stars in ten star clusters, as presented by Deason et al. (2011). Unlike Deason et al. (2011), Barbosa et al. (2022) argue that

$$M_{g_{SDSS,0},BHB}((g-r)_{SDSS,0}) = \frac{0.178}{0.537+(g-r)_{SDSS,0}} \quad (1)$$

results in smaller uncertainties for the absolute magnitude of BHB stars in Sloan $g_0$ when revisiting the approximation by Deason et al. (2011) with more data of BHB stars in star clusters and better constraints on the distances of the star clusters. In both Deason et al. (2011) and Barbosa et al. (2022), the intrinsic spread in the absolute magnitude-colour relation is of 0.1 mag. This spread takes into account the metallicity-dependence of the colour-absolute magnitude relation in case of Deason et al. (2011). In contrast, Barbosa et al. (2022) find that Equation (1) does not have any significant metallicity-dependence.

We adopt Equation (1) to calculate the heliocentric distances of the BHB stars while noting that the output obtained from the Deason et al. (2011) polynomial are consistent with the distances resulting from Equation (1) within their uncertainties. It is important to stress here that not all the stars used here have an SDSS photometric measurement, in particular for the stars selected from the Thomas et al. (2018) catalogue. In that case, the colour in the SDSS photometric system is computed from the PS1 measurement using the Equation (6) of Thomas et al. (2018). To enable the computation of SDSS $g_0 - M_{SDSS,0}$, $g_{SDSS,0}$ were approximated with $g_{PS,0}$. The uncertainties on the distances were computed by adding in quadrature the uncertainties due to the photometric uncertainties in the g and r bands to the intrinsic scatter of the colour-absolute magnitude relation. This results in uncertainties of typically 5 per cent.

## 3. Results

### 3.1. Comparison between observations and N-body simulations of stars in the distant part of the Sagittarius tidal stream

We explore the relationship between heliocentric distance and (corrected) LOS velocity for the BHB stars in our sample relative to predictions of N-body simulations of the formation and evolution of the Sagittarius stream.

#### 3.1.1. Radial velocity structure of blue horizontal branch stars

Figure 6 examines the radial velocity structure (stars and diamonds for results that are velocity-corrected and not, respectively) of the apocentre of the trailing arm and spur feature of the Sagittarius stream as traced by our sample of BHB stars<fnref n="6" />. The constraints on the kinematics of the distant part of the trailing arm of the Sagittarius stream that is traced by our sample of

<footnote n="6">⁶ For the purpose of Galactic Standard of Rest LOS velocity calculation, we adopted the definition of the Galactocentric coordinate system as also used in Vasiliev et al. (2021) that is the standard in ASTROPY (Astropy Collaboration 2013, 2018, 2022). This definition uses the J2000 International Celestial Reference System coordinates (RA = $17^h45^m37.224$, Dec = $-28^\circ56'10.''23$) of the Galactic centre noted in Reid & Brunthaler (2004), a distance = 8.122 kpc between the Sun and Galactic centre as constrained by GRAVITY Collaboration (2018), height = 20.8 pc of the Sun above the Galactic mid-plane, and (12.9, 245.6, 7.78) km s⁻¹ (Drimmel & Poggio 2018) as the 3D Solar velocity relative to the Galactic centre.</footnote>
<figure><img src="image_9.png" /><figcaption>**Fig. 6.** Comparison of the observational data on the LOS velocities in the Galactic Standard of Rest of the BHB stars (diamonds and stars) and predictions from the Dierickx & Loeb (2017) (light grey points) and Vasiliev et al. (2021) without (dark grey points) and with LMC (orange points) models in the same angular range across the stream. Left panel: big diamond markers indicate LOS velocities in Galactic Standard of Rest of BHB stars estimated with RVSPECFIT from the FORS2 spectra. Right panel: stars represent the same velocities with an applied shift as reported in Figs. 3 and C.1. The plot also shows the position of globular cluster NGC 2419 (filled circle) in this projection of phase space. The LOS velocity of NGC 2419 comes from Vasiliev & Baumgardt (2021) also assuming the right ascension and declination provided in the Vasiliev & Baumgardt (2021) catalogue, while we use the heliocentric distance of NGC 2419 from Baumgardt & Vasiliev (2021).</figcaption></figure>

BHB stars are compared to the present-day data of stellar particles in three simulations of the infall of the Sagittarius dwarf galaxy.

The three simulations we consider are by Dierickx & Loeb (2017, predictions in Fig. 6 are shown with grey points)<fnref n="7" /> and Vasiliev et al. (2021, one takes into account the effects of the LMC and the other does not). Data of the stellar particles in the N-body simulations by Vasiliev et al. (2021) are represented by orange and dark grey points in Fig. 6 for the case of the simulation without and with the LMC, respectively<fnref n="8" />. The main differences between these simulations are as follows:

- no LMC perturbations in Dierickx & Loeb (2017) simulation;

- Hernquist (1990) model for gravitational potential of Galactic bulge with total mass of $1.25 \times 10^{10} M_{\odot}$ in Dierickx & Loeb (2017) versus an exponentially truncated, spheroidal power law model with total mass of $1.2 \times 10^{10} M_{\odot}$ in Vasiliev et al. (2021) for the bulge density;

- exponentially declining gravitational potential of disc with total mass of $8.125 \times 10^{10} M_{\odot}$ in Dierickx & Loeb (2017) versus isothermal sheet with exponential radial variation for the disc density and total mass of $5 \times 10^{10} M_{\odot}$ in Vasiliev et al. (2021);

- spherical Hernquist (1990) model for gravitational potential of Galactic halo with total mass of $1.25 \times 10^{12} M_{\odot}$ in Dierickx & Loeb (2017) versus an exponentially truncated, Zhao (1996) model with varying flattening and orientation in Vasiliev et al. (2021) for halo density;

- adoption by Dierickx & Loeb (2017) of the same models assumed to describe the gravitational potential of the Milky Way for the Sagittarius progenitor where bulge, disc, and halo have total masses of $5.2 \times 10^8 M_{\odot}$, $7.8 \times 10^8 M_{\odot}$, and $1.3 \times 10^{10} M_{\odot}$, while Vasiliev et al. (2021) assume the King

model for the stellar density of the Sagittarius progenitor and a more extended, spherical, cored DM halo of mass of $3.6 \times 10^9 M_{\odot}$ that is gradually tidally disrupted with a mass loss described by a piecewise linear function.

In Fig. 6 there is a clear trend of increasing Galactic Standard of Rest LOS velocities beyond heliocentric distance $D_{\odot} \approx 120$ kpc for the data of the final (present-day) snapshot of the Dierickx & Loeb (2017) and Vasiliev et al. (2021) N-body simulations of the formation and evolution of the Sagittarius stream. To a lesser extent we see a similar trend in the BHB stars (both the estimates of the LOS velocity without and with shifts to correct the systematics outlined in Sect. 2.6 in the left and right panel, respectively). The correlation between heliocentric distance and Galactic Standard of Rest LOS velocity for $D_{\odot} \gtrsim 120$ kpc is interesting because these are the stars expected in the spur feature of the trailing arm of the Sagittarius stream and the fact that the BHB stars trace such a feature qualitatively confirms its existence. The models shown in Fig. 6 and also Fardal et al. (2019)<fnref n="9" /> have significant different predictions for the velocity structure of this feature.

In both panels of Fig. 6, the trend for the spur feature as traced by the BHB stars follows the models from Vasiliev et al. (2021) more closely. Depending on whether the tentative correction is applied (right panel), it follows more the Vasiliev et al. (2021) model with the LMC. However, we also note that in both cases (without and with the correction) the apocentre data points qualitatively fit better with the Vasiliev et al. (2021) model without the LMC.

In Fig. 6, we note that the shifts applied to the LOS velocities of the BHB stars based on the values in Fig. 3 (from the left to the right panel) are not significant compared to the large velocity range that we study. Second-order variations in the LOS velocity measurements such as due to the flexure therefore reasonably do not affect our conclusions. Over half of the BHB stars surveyed in our sample could be part of the apocentre of the trailing arm which is visible in Fig. 6 below 120 kpc as the arc characteristic of the apocentres where stars turn around and have smaller and smaller absolute LOS velocities. Interestingly, the position of the apocentre of the trailing arm according to our

<footnote n="7">⁷ The data was gathered from mdierick.github.io/Sgr_stars_kin_data.csv. The LOS velocities of the particles in their simulation are already in the Galactic Standard of Rest with a slightly different definition of the Galactocentric coordinate system.</footnote>

<footnote n="8">⁸ We use the same procedure for the computation of the LOS velocities in the Galactic Standard of Rest for the stellar particles in the Vasiliev et al. (2021) model as for the BHB stars.</footnote>

<footnote n="9">⁹ The data for this model is not publicly available.</footnote>
data was observed to be systematically off of any the models presented in Fig. 6. In summary, the constrained LOS velocity structure of the apocentre of the trailing arm and the spur beyond it of the Sagittarius stream with our sample of BHB stars is informative to constrain current models of the Sagittarius stream and also in turn the Galactic gravitational potential at these distances.

We note that Li et al. (2023a) present a set of (candidate) M-type giant stars that also reaches the spur region of the distant Sagittarius stream trailing arm but with very limited accuracy in their distances (the uncertainties of the photometric, heliocentric distances $\Delta D_{\odot} \ge 30$ kpc at $D_{\odot} \ge 70$ kpc, hence 3–7 times larger than we find for the BHB stars, see also Appendix D). Our sample traces this region combining accurate distances and velocity space for the first time. In the next two sections we discuss more quantitatively how we can constrain models of the formation and evolution of the Sagittarius stream and the shape of the Galactic gravitational potential at the radii that our stars cover.

### 3.1.2. Distance rescaling of Sagittarius apocentre and spur in Vasiliev et al. (2021) simulations

As discussed above, we find there is a difference in the mapped velocity structure of the BHB stars in the most distant part of the trailing arm of the Sagittarius stream and predictions from the Dierickx & Loeb (2017) and Vasiliev et al. (2021) models (see Fig. 6). Of the two models, the latter resembles the trend expressed by our data in this study more. This model is also a bit closer to the observed position of the apocentre of the trailing arm of the Sagittarius stream in the 2D projection of heliocentric distance and LOS velocity. Nonetheless, this apocentre turn-around seems somewhat offset in distance in the model predictions compared to our data points. A possible explanation for this might be that the enclosed Galactic mass at the apocentre of the trailing arm (here 100 kpc) found by Vasiliev et al. (2021) of $(5.6 \pm 0.4) \times 10^{11} M_{\odot}$ is lower. In consideration of the fact that the Vasiliev et al. (2021) model accurately reproduces the majority of the Sagittarius stream at closer Galactic distances, it is also proposed that the assumed outer slope $\beta$ of the halo density profile by Vasiliev et al. (2021, $\rho \propto (R/\text{rscale})^{-\gamma} (1+(R/\text{rscale})^{\alpha})^{(\gamma-\beta)/\alpha} \exp(-(R/200 \text{kpc})^2)$ and $R \equiv (pq)^{1/3} \sqrt{X^2+(Y/p)^2+(Z/q)^2}$) is not an accurate description of the true potential at these (relatively unconstrained) distances.

To investigate how $\beta$ and in turn the enclosed mass of the MW within 100 kpc has to be adapted to match our constraints, we rescale across heliocentric distance (multiplicative factor $f$ for the heliocentric distance) the data of the stellar particles associated to the trailing arm region of the Sagittarius stream in the presented N-body simulations by Vasiliev et al. (2021) without and with the LMC (cf. with Fig. 6). Average tracks were generated for all data of the Vasiliev et al. (2021) models in Fig. 6 (stellar particles at $-230^{\circ} < \Lambda_{Sgr} < -160^{\circ}$ in the final, present-day snapshots of the N-body simulations). These tracks are shown in Fig. 7, where the top panel points were obtained by grouping the 500 closest data points together in heliocentric distance and velocity space using K-means clustering as implemented in SCIKIT-LEARN (Pedregosa et al. 2011) and finding the standard deviations that describe each group/cluster of data points.

Chi-square

$$ \chi^2 \approx -0.5 \sum_{k=1}^{K} \frac{(f \mu_{D_{\odot},k} - D_{\odot,k})^2}{(f \sigma_{D_{\odot},k})^2 + \Delta D_{\odot,k}^2} + \frac{(\mu_{V_{h,k}} - V_{h,k})^2}{\sigma_{V_{h,k}}^2 + \Delta V_{h,k}^2} \quad (2) $$

<figure><img src="image_10.png" /><figcaption>Fig. 7. Rescaling of the average track of stellar particles in heliocentric distance and LOS velocity in final (present-day) snapshots of N-body simulations of the Sagittarius stream by Vasiliev et al. (2021) without and with perturbations by the LMC across heliocentric distance. (a) Copy of the right panel of Fig. 6 except replacing Galactic Standard of Rest by heliocentric LOS velocities, removing the data of the particles in the Dierickx & Loeb (2017) model, and adding average tracks of the data of the stellar particles in the Vasiliev et al. (2021) N-body simulations. (b) heliocentric 2D phase space data of stellar particles of N-body simulations presented in Vasiliev et al. (2021) after rescaling to match it to the shown constraints from the BHB stars.</figcaption></figure>

tests were used to compare the closest nodes $(f \mu_{D_{\odot},k} \pm f \sigma_{D_{\odot},k}, \mu_{V_{h,k}} \pm \sigma_{V_{h,k}})$ of the average tracks to four considered cases of some representation of the K data points $(D_{\odot,k} \pm \Delta D_{\odot,k}, V_{h,k} \pm \Delta V_{h,k})$ in Figs. 6 and 7 of each star excluding outlier candBHB4$^{10}$. These four cases are

- (B22, Heliocentric) with heliocentric distance of BHB stars estimated from Equation (1) and LOS velocities of BHB stars (data represented in the left panel of Fig. 6).

- (B22, Corrected) with heliocentric distance of BHB stars estimated from Equation (1) and corrected LOS velocities of BHB stars (data represented in the right panel of Fig. 6).

- (D11, Heliocentric) heliocentric distance of BHB stars estimated from Deason et al. (2011) relation and LOS velocities of BHB stars.

- (D11, Corrected) with heliocentric distance of BHB stars estimated from Deason et al. (2011) relation and corrected LOS velocities of BHB stars.

Subsequently, we sample the natural logarithm of the posterior of the rescaling factor $f$ given the data that we calculated from the natural logarithm of a Gaussian likelihood (equivalent to approximately $\chi^2$) of the difference between the data of the BHB

<footnote>$^{10}$ No significant difference were found between results from these tests with and without the data of candBHB4. Its classification is uncertain (see Table 3).</footnote>
<table><caption>**Table 4.** Differences between the data of the velocity structure of the BHB stars at the apocentre and spur feature of the trailing arm of the Sagittarius stream to which the data of the final (present-day) snapshots of the N-body simulations of the formation and evolution of the Sagittarius stream in the presence or absence of the LMC by Vasiliev et al. (2021) were scaled along heliocentric distance to match it.</caption><thead><tr><th>Label<br>(1)</th><th>Dist.<br>(2)</th><th>LOS vel.<br>(3)</th><th>LMC<br>(4)</th><th>f<br>(5)</th><th>χ<sup>2</sup><br>(6)</th></tr></thead><tbody><tr><td>A</td><td>D11</td><td>Heliocen.</td><td>No</td><td>1.04<sup>+0.02</sup><sub>-0.04</sub></td><td>-16.3<sup>+0.3</sup><sub>-0.9</sub></td></tr><tr><td>B</td><td>B22</td><td>Heliocen.</td><td>No</td><td>1.05<sup>+0.02</sup><sub>-0.03</sub></td><td>-16.4<sup>+0.2</sup><sub>-0.9</sub></td></tr><tr><td>C</td><td>D11</td><td>Corrected</td><td>No</td><td>1.03<sup>+0.02</sup><sub>-0.02</sub></td><td>-33.7<sup>+0.2</sup><sub>-1.1</sub></td></tr><tr><td>D</td><td>B22</td><td>Corrected</td><td>No</td><td>1.05<sup>+0.02</sup><sub>-0.02</sub></td><td>-33.2<sup>+0.2</sup><sub>-0.7</sub></td></tr><tr><td>E</td><td>D11</td><td>Heliocen.</td><td>Yes</td><td>1.12<sup>+0.02</sup><sub>-0.02</sub></td><td>-8.4<sup>+0.3</sup><sub>-0.9</sub></td></tr><tr><td>F</td><td>B22</td><td>Heliocen.</td><td>Yes</td><td>1.15<sup>+0.02</sup><sub>-0.02</sub></td><td>-10.0<sup>+0.2</sup><sub>-0.5</sub></td></tr><tr><td>G</td><td>D11</td><td>Corrected</td><td>Yes</td><td>1.10<sup>+0.02</sup><sub>-0.02</sub></td><td>-12.6<sup>+0.5</sup><sub>-0.5</sub></td></tr><tr><td>H</td><td>B22</td><td>Corrected</td><td>Yes</td><td>1.11<sup>+0.02</sup><sub>-0.02</sub></td><td>-11.7<sup>+0.2</sup><sub>-1.1</sub></td></tr></tbody></table>

**Notes.** The second column is the abbreviation for the equation that we used to calculate the absolute magnitude in Sloan *g*<sub>SDSS</sub> of the sample BHB stars, and as a result the distance modulus and heliocentric distance of the stars. D11 refers to Deason et al. (2011); B22 refers to Barbosa et al. (2022). Column (3) lists LOS velocity of each BHB star in the set that we study we consider for the rescaling: ‘Heliocen.’ refers to the heliocentric LOS velocity without a shift to correct for systematics as estimated from the velocity standard stars (see Fig. 3 or Fig. C.1 and left panel of Fig. 6); ‘Corrected’ refers to the corrected (heliocentric) LOS velocities (see right panel of Fig. 6). Column (4) indicates whether this is the data from the N-body simulation by Vasiliev et al. (2021) with (‘Yes’) or without (‘No’) perturbations by the LMC. Column (6) gives the approximate χ² of the analysis as calculated from Equation (2) (more precisely the natural logarithm of the assumed Gaussian likelihood).

stars and the average tracks to find the optimal value of *f*. The uncertainties of the average tracks and observational data were added in quadrature but the resulting value of *f* is dominated by the uncertainties of the observational measurements, rather than by the uncertainties of the average track fitting, meaning that the exact fitting procedure (e.g. the number of nodes chosen) does not affect the results significantly. We adopt these uncertainties for all values derived. We adopt a uniform prior for *f* over 0.1 and 10.0. This sampling across Markov chains with a Monte Carlo algorithm and optimisation was achieved with EMCEE (Foreman-Mackey et al. 2013) using ten walkers and 1000 steps. Table 4 lists the results of the rescaling (median and 16 and 84 quantiles of sampled posterior distribution of *f*). We note that none of the rescaling experiments within the group of cases including the LMC influence (‘E’, ‘F’, ‘G’, and ‘H’) as well as within the group of cases without the LMC influence (‘A’, ‘B’, ‘C’, and ‘D’) resulted in differences that were statistically significant. This illustrates again that – considering the main results presented in this paper – the possible flexure corrections of the radial velocities are not that influential and, moreover, that the choice of distance calibration for the BHB stars does not alter the results significantly. This minimal effect of the flexure corrections and choice of distance calibration is fully expected as the either choice stays within the uncertainties of these parameters. On the other hand, as maybe expected, the difference for *f* between the two groups (with and without the LMC perturbation) is statistically significant. For the remainder of this work, case ‘D’ (B22, Corrected, No) and ‘H’ (B22, Corrected, Yes) are labeled as our fiducial cases with and without the LMC.

Panel b of Fig. 7 compares the rescaling of the 2D heliocentric phase space (distance, LOS velocity) data of the stellar particles at $-230^{\circ} < \Lambda_{Sgr} < -160^{\circ}$ along the Sagittarius stream in the present-day snapshots of the N-body simulations without ($f = 1.05^{+0.02}_{-0.02}$ for fiducial case) and with the LMC ($f = 1.11^{+0.02}_{-0.02}$ for fiducial case) as performed by Vasiliev et al. (2021) for the fiducial cases. The results indicate that the rescaled Vasiliev et al. (2021) model with the LMC fits the data better than the no LMC model quantitatively, given the $\chi^2$ values from the analysis (see Equation (2)) in the last column of Table 4. This supports a number of previous works that demonstrate the effects of the LMC on the MW and the over-densities in it (e.g. Erkal et al. 2019; Garavito-Camargo et al. 2019; Shipp et al. 2019; Petersen & Peñarrubia 2021; Vasiliev et al. 2021; Garavito-Camargo et al. 2021; Lilleengen et al. 2023; Koposov et al. 2023; Chandra et al. 2024; Byström et al. 2024).

In addition, the shift in the distance modulus as a result of the rescaling is $5 \log_{10} f = 0.23 \pm 0.04$ mag. This shift is significantly larger than the spread observed for the colour-$M_{g_{SDSS,0}}$ used to estimate the heliocentric distances to the BHB stars.

## 3.2. Milky Way mass profile

Given the similarities between the rescaled Vasiliev et al. (2021) model with the LMC and the observed Sagittarius stream, it is possible to use this model to re-evaluate the MW’s mass within the Sagittarius apocentre radius ($\approx 100$ kpc). In this section, we present a first analysis in this direction to provide an estimate of the (magnitude of) change in mass implied, but we note that a full and detailed re-modelling is beyond the scope of this work.

The rescaling factor derived in Section 3.1.2 is computed such that the heliocentric LOS velocities in the model match the observed, heliocentric velocities at the rescaled Galactocentric positions<fnref n="11" />. In this analysis, we exclusively utilise the LOS velocity. Note that proper motions are not considered because they are not available for all stars.

In our approach, we assume that the Sagittarius stream model from Vasiliev et al. (2021) provides a good representation of the total energy distribution at every position along the stream. However, as shown in Fig. 7, we have seen that the model needs to be stretched to match the distance and velocity distribution of the data. Therefore, we stretch the Galactocentric coordinates of the Sagittarius stream model of Vasiliev et al. (2021) (with the LMC) according to the rescaling factor found in Section 3.1.2, but leave the 3D velocities of the particles unchanged. Adjusting the particles in the Vasiliev et al. (2021) simulation according to this rescaling factor provides direct access to the gravitational potential ($\Phi$) of each star given that $E \equiv \Phi + 0.5 V^2$, where E is the total energy and V the Galactocentric velocity of the star.

Consequently, we adjust the MW DM halo profile such that the gravitational potential with the new profile, $\Phi_{MW,\text{new}}(R)$, of a star at the rescaled Galactocentric position $R'$ matches the potential of the Vasiliev et al. (2021) model, $\Phi_{MW,\text{Vasiliev+2021}}$, at Galactocentric position $R$ before rescaling.

Moreover, given that the original (non-rescaled) Vasiliev et al. (2021) model accurately reproduces most of the Sagittarius stream at shorter Galactic distances, we only adjust the outer slope of the MW DM halo. Specifically, we assume that the triaxiality, twisting, and baryonic distribution remain unchanged

<footnote n="11">¹¹ While the rescaling results presented in Section 3.1.2 are across heliocentric distances, we use here the effect of the rescaling of the heliocentric distances on the Galactocentric positions and furthermore note that rescaling in Galactocentric distances instead yields consistent results compared to the entries in Table 4.</footnote>
compared to the MW model proposed by Vasiliev et al. (2021)<fnref n="12" />. In particular, we note that we thereby assume that the gravitational perturbation of the LMC remains unchanged at the rescaled position. This assumption is not entirely correct, as for a same fixed mass of the LMC as found by Vasiliev et al. (2021), its perturbation on the apocentre of the Sagittarius stream in the rescaled model should be less important than at the original positions of the Vasiliev et al. (2021) model, given that the rescaled position are further away from the LMC than the original position. At the same time, the mass ratio of the LMC and MW will be more comparable with a lower mass MW. Properly accounting for these perturbations would require rerunning the simulation of Vasiliev et al. (2021), as the LMC directly exerts a gravitational force on the Sagittarius stream but also modifies the shape of the DM halo (Garavito-Camargo et al. 2019; Petersen & Peñarrubia 2021). However, we stress this will be a secondary effect due to the relatively small changes in distance compared to the original model.

We impose the boundary condition on the gravitational potential that the gravitational potential at the Galactocentric location of the Sagittarius dwarf spheroidal galaxy remnant, $R_{Sgr} = (17.9, 2.6, -6.6)$ kpc, should be the same between the original Vasiliev et al. (2021) model and the adjusted DM profile:

$$ \Phi_{\text{MW,new}}(\mathbf{R}_{\text{Sgr}}) = \Phi_{\text{MW,Vasiliev+2021}}(\mathbf{R}_{\text{Sgr}}). $$

Thus, revising the enclosed Galactic mass within 100 kpc and by extension exploring the changes in the outer slope $\beta$ of the MW DM halo involves minimising the following least-squares function:

$$ \chi^2 = \sum_{i=1}^{N} [\Phi_{\text{MW,new}}(\mathbf{R}'_i) - \Phi_{\text{MW,Vasiliev+2021}}(\mathbf{R}_i) - \text{norm}]^2. \quad (3) $$

Here $\mathbf{R}'_i$ represents the rescaled Galactocentric positions along the stream and norm accounts for normalisation adjustments. The subscript $i$ represents the $N$ particles in the distant trailing arm region of the Sagittarius stream ($-230^\circ < \Lambda_{Sgr} < -160^\circ$) from the Vasiliev et al. (2021) simulation. In the region of the Sagittarius stream with the data of the BHB stars, the positional, Galactocentric vector $(X_i, Y_i, Z_i)$ is represented by $\mathbf{R}_i$ for the $i$th particle. $\Phi_{\text{MW,new}}$ is calculated using the stretched, Galactocentric $X'_i, Y'_i$, and $Z'_i$ positions found from the rescaling described in Sect. 3.1 to yield the potential in $\mathbf{R}'_i = (X'_i, Y'_i, Z'_i)$ for the $i$th particle. All potentials are computed with the Action-based Galaxy Modelling Architecture code (AGAMA, Vasiliev 2018, 2019a) version 1.0.

We employ a Monte Carlo technique to sample from the unknown distribution of the enclosed MW mass within 100 kpc and $\beta$ by repeating the minimisation of Equation (3) to find the optimal value of the mass and $\beta$ with 1000 realisations of $f$ and thus $\mathbf{R}'_i$ and the other parameters of the halo density profile model by Vasiliev et al. (2021). In the Monte Carlo method, for each realisation we vary the rescaling factor $f = 1.11 \pm 0.02$, scale radius rscale, inner slope $\gamma$, and transition steepness $\alpha$ of the halo density profile model by Vasiliev et al. (2021) (Eugene Vasiliev, private communication) and the other parameters of the model by Vasiliev et al. (2021) that quantify the triaxiality and twisting of the halo (see Figure 14 in Vasiliev et al. 2021) within their uncertainties. We use split normal distributions to approximate

<figure><img src="image_11.png" /><figcaption>Fig. 8. Estimates of enclosed mass of Milky Way at different radii together with our constraints with the (blue) star symbol at 100 kpc. As in the diagram in Bobylev & Baykova (2023), which we adapted here, with the exceptions of 32 and 33 which are new estimates by Ibata et al. (2024) and are highlighted here, the annotated numbers for each data point correspond to 1–Küpper et al. (2015), 2–Malhan & Ibata (2019), 3–Prudil et al. (2022), 4–Posti & Helmi (2019), 5–Williams et al. (2017), 6–Ablimit & Zhao (2017), 7–Williams & Evans (2015), 8–Xue et al. (2008), 9–Gnedin et al. (2010), 10–Gibbons et al. (2014), 11–Vasiliev et al. (2021), 12–Shen et al. (2022), 13–Correa Magnus & Vasiliev (2022), 14–Eadie & Jurić (2019), 15–McMillan (2017), 16–Vasiliev (2019b), 17–Battaglia et al. (2005), 18–Eadie et al. (2017), 19–Eadie & Harris (2016), 20–Deason et al. (2012), 21–Ablimit et al. (2020), 22–Bird et al. (2022) (from data of K-type giants), 23–Zhou et al. (2023), 24–Bhattacharjee et al. (2014), 25–Wang et al. (2022), 26–Bird et al. (2022) (from data of BHB stars), 27–Sun et al. (2023), 28–Bajkova & Bobylev (2016), 29–Huang et al. (2016), 30–Eadie et al. (2015), 31–Patel et al. (2018), 32–Ibata et al. (2024) (within 50 kpc), and 33–Ibata et al. (2024) (estimate of virial mass within 200 kpc).</figcaption></figure>

the distributions of the parameters of the halo density profile model by Vasiliev et al. (2021). We do not take into account here that some of the parameters of the Vasiliev et al. (2021) model are correlated with $\beta$, but treat them as independent. The methodology employed here is relatively straightforward and, it must be acknowledged, somewhat simplistic. It is assumed that the energy distribution for the Sagittarius progenitor provided by the Vasiliev et al. (2021) model is reliable, given that it already closely matches the majority of the stream. We further note here that the conserved energy distribution combined with altered distances will lead to a corresponding relative change in angular momentum compared to the original model (see Figure 11 of Vasiliev et al. 2021). However, it is our contention that it can provide a satisfactory initial approximation of the local potential before undertaking further, more comprehensive simulations in future work.

The results of the sampling show that $\beta = 2.47_{-0.21}^{+0.26}$ can lead to a position of the apocentre of the trailing arm of the Sagittarius stream in the Vasiliev et al. (2021) that is in line with our constraints based on the BHB stars. This value of $\beta$ is slightly higher, but within the uncertainties of the $\beta$ value reported by Vasiliev et al. (2021, $\beta = 2.43_{-0.13}^{+0.16}$), Eugene Vasiliev, private communication). Consequently, the enclosed mass at the apocentre of the trailing arm (here ~100 kpc) by Vasiliev et al. (2021) is lowered to $(5.3 \pm 0.4) \times 10^{11} M_\odot$, again consistent within the uncertainties of Vasiliev et al. (2021) who report $M(< 100 \text{ kpc}) = (5.6 \pm 0.4) \times 10^{11} M_\odot$.

Figure 8 shows our new estimate of the enclosed mass of the Milky Way within 100 kpc in comparison to other estimates of

<footnote n="12">¹² Our adjustment of the MW DM profile incorporates uncertainties in these parameters.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Proper motion of the identified BHB stars (stars) in comparison to the simulation data by Dierickx & Loeb (2017) (point markers) and Vasiliev et al. (2021) (square and plus markers showing rescaled proper motions by $1/f$, where $f$ is the same factor used to rescale the heliocentric distances of the simulation data, i.e. $1.11 \pm 0.02$; we note that the uncertainties in the rescaled proper motions resulting from the uncertainty in $f$ are smaller than the symbol size) and the globular cluster NGC 2419 (filled circle). We show the predicted motion along right ascension and declination in the simulations across the same angular range of the stream as the BHB stars cover. The proper motion estimates of NGC 2419 correspond to the values provided by Vasiliev & Baumgardt (2021).</figcaption></figure>

the enclosed mass at various radii between 10 kpc and 270 kpc. Our results are in line with the general trend of the literature.

### 3.3. NGC 2419

Finally, as outlined in the introduction, we use our new dataset to investigate more closely the proposed connection between the trailing arm of Sagittarius stream and the globular cluster NGC 2419 (e.g. Irwin 1999; Newberg et al. 2003; Ruhland et al. 2011; Belokurov et al. 2014; Sohn et al. 2018; Massari et al. 2019; Antoja et al. 2020; Bellazzini et al. 2020; Peñarrubia & Petersen 2021; Davies et al. 2024; Chen & Gnedin 2024; Rostami Shirazi et al. 2024). The (mean) position of NGC 2419 in Fig. 6 was found to be consistent with the approaching group of BHB stars at the apocentre of the trailing arm.

Additionally, a subsample of BHB stars with estimates of proper motion in Gaia DR3 was prepared. Here, we adapt the procedure suggested by Fabricius et al. (2021) and deselect spurious astrometric solutions in Gaia DR3 with $IPD\_GOF\_HARMONIC\_AMPLITUDE \ge 0.1$. This excludes candBHB8 and candBHB16. Fig. 9 shows that the remaining targets with Gaia DR3 proper motions across right ascension and declination (Gaia Collaboration 2023) tentatively agree with data for NGC 2419 by Vasiliev & Baumgardt (2021) as well as with model predictions by Dierickx & Loeb (2017) and Vasiliev et al. (2021). CandBHB9, candBHB10, candBHB11, candBHB13, candBHB17, candBHB19, and candBHB21 could therefore be linked to the apocentre of the trailing arm of the Sagittarius stream and its extension, the spur feature, across phase space.

## 4. Conclusion

The objectives of this study are twofold: to evaluate the effectiveness of selecting candidate BHB stars in the outer halo beyond

~80 kpc and to constrain the kinematics of the apocentre and spur feature of the sparsely studied outer regions of the trailing arm of the Sagittarius stream. Our ESO/VLFT/FORS2 pilot programme confirmed that 20 out of the 25 candidate BHB stars chosen from the samples by Thomas et al. (2018) and Starkenburg et al. (2019) are true BHB stars based on the spectroscopically estimated effective temperatures and surface gravities from fitting with RVSPECFIT (Fig. 4). As a consequence, we can confidently push the photometric selection techniques developed in Thomas et al. (2018) and Starkenburg et al. (2019) to stars fainter than the limits of Gaia (G ≈ 21 mag, Gaia Collaboration 2023). The present study has been one of the first attempts to examine the performance of photometric selections of candidate BHB stars at these distances with a dedicated spectroscopic follow-up program.

Despite its exploratory nature, this study offers also some insight into the behaviour on the VLT/FORS2 instrument and the 600B+22 grism we used when observing stars at high air mass (above 1.5) where second-order effects in the mapping from pixel to wavelength space across the detector are expected due to telescope flexure under gravity. Although the current study is based on a small sample of candidate BHB and velocity standard stars to assess and quantify the above-mentioned variation(s), the findings from adopting the same methodology to estimate heliocentric LOS velocities of the velocity standard stars as employed for the target stars suggest systematic offsets from reference values of the order of at least $10–25 \text{ km s}^{-1}$. These expected and corrected systematics in the LOS velocity measurements of the target stars in our observational set-up made no significant difference to the overall conclusions that we draw from the available phase space data of the BHB stars (Sects. 2.6 and 6).

The FORS2 spectra of the confirmed BHB stars provide further observational constraints on the current Milky Way and Sagittarius stream models (e.g. Dierickx & Loeb 2017; Fardal et al. 2019; Vasiliev et al. 2021; Hainje et al. 2025), using the selected data on the spur feature of the Sagittarius debris. The velocity structure of the apocentre of the trailing arm of the Sagittarius stream, at heliocentric distances of $D_{\odot} \approx 100 \text{ kpc}$, as traced by the identified BHB stars, is offset from the model predictions by Dierickx & Loeb (2017) and Vasiliev et al. (2021, Fig. 6). The investigation of heliocentric distance and LOS velocity in the N-body simulation of the disruption of the Sagittarius dwarf galaxy in the presence of the LMC by Vasiliev et al. (2021) in comparison to the data of the BHB stars has shown that the observations can be described well by the model predictions by Vasiliev et al. (2021), including the spur feature, if the heliocentric distances of the simulation data are rescaled by a factor $f = 1.11_{-0.02}^{+0.02}$. This would correspond to a change in the outer slope of $\beta = 2.47_{-0.21}^{+0.26}$ (versus $2.43_{-0.13}^{+0.16}$) and in turn leads to a refined estimate of the enclosed Milky Way mass within the apocentre of the trailing arm that is slightly lower but consistent within the uncertainty compared to the estimate by Vasiliev et al. (2021) at $(5.3 \pm 0.4) \times 10^{11} \text{ M}_{\odot}$.

A limitation of this study is that the use of BHB stars as tracers of the metal-poor Milky Way halo does not permit the investigation of metallicities. This is due to the difficulty of constraining the metallicities of hot stars in general through optical spectroscopy, and in particular in the cases of low spectral resolution, due to the blending of metallicity-sensitive lines by strong Balmer lines. Future studies combining different stellar tracers where estimates of metallicity can also be obtained (such as for instance, red giants, or RR Lyrae; e.g. Viswanathan et al. 2024; Medina et al. 2025) will shed more light on the kinematics and
also on the chemistry of distant substructures in the Galactic halo.

The current data highlight the importance of providing a more detailed view of the outer Galactic halo through various standard candle stars that also serve as good kinematic tracers. The findings of this research provide first insights on the kinematics of one of the most distant extensions of the Sagittarius stream currently known, the spur feature. We expect that the insights gained from this study will provide helpful constraints to modellers of the Sagittarius stream and the Galactic potential.

**Acknowledgements.** The authors like to thank the anonymous referee for a careful and constructive report that has certainly helped to improve the manuscript. We want to express our gratitude to Eugene Vasiliev for providing us the uncertainties of the scale radius, inner and outer slope, and transition steepness of the fiducial model of the halo density profile in Vasiliev et al. (2021). MB would also like to thank Sergey E. Koposov for their help in understanding and using some of the new features in RVSPECFIT (Koposov et al. 2011; Koposov 2019). MB and ES acknowledge funding through VIDI grant "Pushing Galactic Archaeology to its limits" (with project number VI.Vidi.193093) which is funded by the Dutch Research Council (NWO). This research has been partially funded from a Spinoza award by NWO (SPI 78-411). GFT and EFA acknowledge support from the Agencia Estatal de Investigación del Ministerio de Ciencia en Innovación (AEI-MICIN) and the European Regional Development Fund (ERDF) under grant number PID2020-118778GB-100/10.13039/501100011033 and the AEI under grant number CEX2019-000920-S. EFA acknowledges support from HORIZON TMA MSCA Postdoctoral Fellowships Project TEMPOS, number 101066193, call HORIZON-MSCA-2021-PF-01, by the European Research Executive Agency. AB acknowledges the Edinburgh Doctoral College Scholarship that funded her work. GEM acknowledges financial support from Natural Sciences and Engineering Research Council of Canada (NSERC) through grant RGPIN-2022-04794, and from an Arts & Sciences Postdoctoral Fellowship at the University of Toronto. AV gratefully acknowledges support from the Canadian Institute for Theoretical Astrophysics (CITA) through a CITA National Fellowship. This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI: 10.26093/cds/vizier). The original description of the VizieR service was published in 2000, A&AS 143, 23. This work is based on data obtained as part of the Canada-France Imaging Survey, a CFHT large programme of the National Research Council of Canada and the French Centre National de la Recherche Scientifique. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA Saclay, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers (INSU) of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project 540 (The Early Milky Way). We additionally made use of ADJUSTTEXT, version 1.2.0 (Flyamer et al. 2024), Action-based Galaxy Modelling Architecture code (AGAMA, Vasiliev 2018, 2019a), ASTROQUERY (Ginsburg et al. 2019), ASTROPY (Astropy Collaboration 2013, 2018, 2022), DUSTMAPS (Green 2018), EMCEE (Foreman-Mackey et al. 2013), ESOREFLEX (Freudling et al. 2013), GALA (Price-Whelan 2017), Image Reduction and Analysis Facility (IRAF, Tody 1986, 1993; National Optical Astronomy Observatories 1999), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), RVSPECFIT (Koposov et al. 2011; Koposov 2019), SCIKIT-LEARN (Pedregosa et al. 2011), SCIPY (Virtanen et al. 2020), SEABORN (Waskom 2021), and VAEX (Breddels & Veljanoski 2018). We would like to acknowledge the use of DeepL Write, an artificial intelligence writing assistant developed by DeepL, which was employed to enhance the readability and style of this article. Author contributions: Conceptualisation, writing and proofreading ESO observing proposals: ES, GFT, NFM, AH, EFA, RI, PJ, GK, TM, RSJ. MB and ES contributed in ESO/VLT/FORS2 data reduction to create the spectral dataset which is employed. Methodology: MB, ES, GFT. Software: MB, ES, GFT. Formal Analysis: MB, ES, GFT. Writing – Original Draft: MB, ES, GFT. All authors helped in finalising the paper draft. Visualisation: MB, ES, GFT, AH, AB, TdB, GK, TM, GEM, SR, RSJ, FS, AV.

## References

* Abdurro'uf, Accetta, K., Aerts, C., et al. 2022, ApJS, 259, 35

* Ablimit, I., & Zhao, G. 2017, ApJ, 846, 10

* Ablimit, I., Zhao, G., Flynn, C., & Bird, S. A. 2020, ApJ, 895, L12

* Almeida, A., Anderson, S. F., Argudo-Fernández, M., et al. 2023, ApJS, 267, 44

* Amarante, J. A. S., Koposov, S. E., & Laporte, C. F. P. 2024, A&A, 690, A166

* Anderson, J. 2020, Very Large Telescope, Paranal Science Operations, FORS2 User Manual, VLT-MAN-ESO-13100-1543

* Antoja, T., Ramos, P., Mateu, C., et al. 2020, A&A, 635, L3

* Astropy Collaboration (Robitaille, T. P., et al.) 2013, A&A, 558, A33

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2018, AJ, 156, 123

* Astropy Collaboration (Price-Whelan, A. M., et al.) 2022, ApJ, 935, 167

* Bailer-Jones, C. A. L. 2017, Practical Bayesian Inference, 142

* Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147

* Bajkova, A. T., & Bobylev, V. V. 2016, Astron. Lett., 42, 567

* Barbosa, F. O., Santucci, R. M., Rossi, S., et al. 2022, ApJ, 940, 30

* Battaglia, G., Helmi, A., Morrison, H., et al. 2005, MNRAS, 364, 433

* Baumgardt, H., & Vasiliev, E. 2021, MNRAS, 505, 5957

* Bellazzini, M., Ibata, R., Malhan, K., et al. 2020, A&A, 636, A107

* Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137

* Belokurov, V., Koposov, S. E., Evans, N. W., et al. 2014, MNRAS, 437, 116

* Bhattacharjee, P., Chaudhury, S., & Kundu, S. 2014, ApJ, 785, 63

* Bird, S. A., Xue, X.-X., Liu, C., et al. 2022, MNRAS, 516, 731

* Bobylev, V. V., & Baykova, A. T. 2023, Astron. Rep., 67, 812

* Breddels, M. A., & Veljanoski, J. 2018, A&A, 618, A13

* Byström, A., Koposov, S. E., Lilleengen, S., et al. 2024, MNRAS, submitted [arXiv:2410.09149]

* Caffau, E., Bonifacio, P., Sbordone, L., et al. 2020, MNRAS, 493, 4677

* Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv e-prints [arXiv:1612.05560]

* Chandra, V., Naidu, R. P., Conroy, C., et al. 2024, ApJ, accepted [arXiv:2406.01676]

* Chen, Y., & Gnedin, O. Y. 2024, Open J. Astrophys., 7, 23

* Clementini, G., Ripepi, V., Garofalo, A., et al. 2023, A&A, 674, A18

* Clewley, L., Warren, S. J., Hewett, P. C., et al. 2002, MNRAS, 337, 87

* Cooper, A. P., Koposov, S. E., Allende Prieto, C., et al. 2023, ApJ, 947, 37

* Correa Magnus, L., & Vasiliev, E. 2022, MNRAS, 511, 2610

* Creevy, O. L., Sordo, R., Pailler, F., et al. 2023, A&A, 674, A26

* Cropper, M., Katz, D., Sartoretti, P., et al. 2018, A&A, 616, A5

* Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Res. Astron. Astrophys., 12, 1197

* Davies, E. Y., Belokurov, V., Monty, S., & Evans, N. W. 2024, MNRAS, 529, L73

* Deason, A. J., Belokurov, V., & Evans, N. W. 2011, MNRAS, 416, 2903

* Deason, A. J., Belokurov, V., Evans, N. W., et al. 2012, MNRAS, 425, 2840

* Deason, A. J., Belokurov, V., Koposov, S. E., & Rockosi, C. M. 2014, ApJ, 787, 30

* Deason, A. J., Belokurov, V., & Koposov, S. E. 2018, ApJ, 852, L18

* DESI Collaboration (Aghamousa, A., et al.) 2016, arXiv e-prints [arXiv:1611.00036]

* Dierickx, M. I. P., & Loeb, A. 2017, ApJ, 836, 92

* Drimmel, R., & Poggio, E. 2018, RNAAS, 2, 210

* Eadie, G. M., & Harris, W. E. 2016, ApJ, 829, 108

* Eadie, G., & Jurić, M. 2019, ApJ, 875, 159

* Eadie, G. M., Harris, W. E., & Widrow, L. M. 2015, ApJ, 806, 54

* Eadie, G. M., Springford, A., & Harris, W. E. 2017, ApJ, 835, 167

* Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, MNRAS, 487, 2685

* Euclid Collaboration (Mellier, Y., et al.) 2025, A&A, 697, A1

* Fabricius, C., Luri, X., Arenou, F., et al. 2021, A&A, 649, A5

* Fantin, N. J., Côté, P., McConnachie, A. W., et al. 2019, ApJ, 887, 148

* Fardal, M. A., van der Marel, R. P., Law, D. R., et al. 2019, MNRAS, 483, 4724

* Flyamer, I., Xue, Z., Colin, et al. 2024, https://doi.org/10.5281/zenodo.12570011

* Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

* Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96

* Fukushima, T., Chiba, M., Homma, D., et al. 2018, PASJ, 70, 69

* Fukushima, T., Chiba, M., Tanaka, M., et al. 2019, PASJ, 71, 72

* Fukushima, T., Chiba, M., Tanaka, M., et al. 2025, PASJ, 77, 178

* Gaia Collaboration (Brown, A. G. A., et al.) 2021, A&A, 649, A1

* Gaia Collaboration (Vallenari, A., et al.) 2023, A&A, 674, A1

* Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2019, ApJ, 884, 51

* Garavito-Camargo, N., Besla, G., Laporte, C. F. P., et al. 2021, ApJ, 919, 109

* Gibbons, S. L. J., Belokurov, V., & Evans, N. W. 2014, MNRAS, 445, 3788

* Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98

* Gnedin, O. Y., Brown, W. R., Geller, M. J., & Kenyon, S. J. 2010, ApJ, 720, L108

* GRAVITY Collaboration (Abuter, R., et al.) 2018, A&A, 615, L15

* Green, G. 2018, J. Open Source Softw., 3, 695

* Gwyn, S., McConnachie, A. W., Cuillardre, J.-C., et al. 2025, AJ, submitted [arXiv:2503.13783]

* Hainje, C., Slone, O., Lisanti, M., & Erkal, D. 2025, arXiv e-prints, [arXiv:2503.15589]

* Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357

* Helmi, A. 2004, ApJ, 610, L97

* Helmi, A., & White, S. D. M. 2001, MNRAS, 323, 529
Hernitschek, N., Sesar, B., Rix, H.-W., et al. 2017, ApJ, 850, 96

Hernquist, L. 1990, ApJ, 356, 359

Huang, Y., Liu, X. W., Yuan, H. B., et al. 2016, MNRAS, 463, 2623

Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90

Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6

Ibata, R. A., & Lewis, G. F. 1998, ApJ, 500, 575

Ibata, R. A., Gilmore, G., & Irwin, M. J. 1994, Nature, 370, 194

Ibata, R. A., Wyse, R. F. G., Gilmore, G., Irwin, M. J., & Suntzeff, N. B. 1997, AJ, 113, 634

Ibata, R., Irwin, M., Lewis, G. F., & Stolte, A. 2001, ApJ, 547, L133

Ibata, R. A., McConnachie, A., Cuillandre, J.-C., et al. 2017, ApJ, 848, 128

Ibata, R., Bellazzini, M., Thomas, G., et al. 2020, ApJ, 891, L19

Ibata, R., Malhan, K., Tenachi, W., et al. 2024, ApJ, 967, 89

Irwin, M. 1999, in IAU Symposium, 192, The Stellar Content of Local Group Galaxies, eds. P. Whitelock, & R. Cannon, 409

Jensen, J., Thomas, G., McConnachie, A. W., et al. 2021, MNRAS, 507, 1923

Ji, A. P., Koposov, S. E., Li, T. S., et al. 2021, ApJ, 921, 32

Johnston, K. V., Law, D. R., & Majewski, S. R. 2005, ApJ, 619, 800

Ju, J., Cui, W., Huo, Z., et al. 2024, ApJS, 270, 11

Koposov, S. E. 2019, Astrophysics Source Code Library [record ascl:1907.013]

Koposov, S. E., Gilmore, G., Walker, M. G., et al. 2011, ApJ, 736, 146

Koposov, S. E., Belokurov, V., Evans, N. W., et al. 2012, ApJ, 750, 80

Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936

Koposov, S. E., Prieto, C. A., Cooper, A. P., et al. 2024, MNRAS, submitted [arXiv:2407.06280]

Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80

Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229

Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008a, AJ, 136, 2022

Lee, Y. S., Beers, T. C., Sivarani, T., et al. 2008b, AJ, 136, 2050

Lenz, D. D., Newberg, J., Rosner, R., Richards, G. T., & Stoughton, C. 1998, ApJS, 119, 121

Li, T. S., Koposov, S. E., Zucker, D. B., et al. 2019, MNRAS, 490, 3508

Li, J., Long, L., Zhong, J., et al. 2023a, ApJS, 266, 4

Li, X.-Y., Huang, Y., Liu, G.-C., Beers, T. C., & Zhang, H.-W. 2023b, ApJ, 944, 88

Lilleengen, S., Petersen, M. S., Erkal, D., et al. 2023, MNRAS, 518, 774

Magnier, E. A., Schlafly, E. F., Finkbeiner, D. P., et al. 2020, ApJS, 251, 6

Majewski, S. R., Skrutskie, M. F., Weinberg, M. D., & Ostheimer, J. C. 2003, ApJ, 599, 1082

Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995

Martin, N. F., Starkenburg, E., Yuan, Z., et al. 2024, A&A, 692, A115

Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4

McMillan, P. J. 2017, MNRAS, 465, 76

Medina, G. E., Li, T. S., Koposov, S. E., et al. 2025, AAS J., submitted [arXiv:2504.02924]

Muraveva, T., Giannetti, A., Clementini, G., Garofalo, A., & Monti, L. 2025, MNRAS, 536, 2749

Naidu, R. P., Conroy, C., Bonaca, A., et al. 2020, ApJ, 901, 48

National Optical Astronomy Observatories 1999, Astrophysics Source Code Library [record ascl:9911.002]

Newberg, H. J., Yanny, B., Grebel, E. K., et al. 2003, ApJ, 596, L191

Oria, P.-A., Ibata, R., Ramos, P., Famaey, B., & Errani, R. 2022, ApJ, 932, L14

Patel, E., Besla, G., Mandel, K., & Sohn, S. T. 2018, ApJ, 857, 78

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, J. Mach. Learn. Res., 12, 2825

Peñarrubia, J., & Petersen, M. S. 2021, MNRAS, 508, L26

Petersen, M. S., & Peñarrubia, J. 2021, Nat. Astron., 5, 251

Posti, L., & Helmi, A. 2019, A&A, 621, A56

Price-Whelan, A. M. 2017, J. Open Source Softw., 2, 388

Prudil, Z., Koch-Hansen, A. J., Lemasle, B., et al. 2022, A&A, 664, A148

Prugniel, P., & Soubiran, C. 2001, A&A, 369, 1048

Ramos, P., Mateu, C., Antoja, T., et al. 2020, A&A, 638, A104

Ramos, P., Antoja, T., Yuan, Z., et al. 2022, A&A, 666, A64

Read, J. I., & Erkal, D. 2019, MNRAS, 487, 5799

Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2023, A&A, 674, A29

Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872

Rostami Shirazi, A., Khalaj, P., & Haghi, H. 2024, MNRAS, 531, 2563

Ruhland, C., Bell, E. F., Rix, H.-W., & Xue, X.-X. 2011, ApJ, 731, 119

Schlafly, E. F., & Finkbeiner, D. P. 2011, ApJ, 737, 103

Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Sesar, B., Hernitschek, N., Dierickx, M. I. P., Fardal, M. A., & Rix, H.-W. 2017a, ApJ, 844, L4

Sesar, B., Hernitschek, N., Mitrović, S., et al. 2017b, AJ, 153, 204

Shen, J., Eadie, G. M., Murray, N., et al. 2022, ApJ, 925, 1

Shipp, N., Li, T. S., Pace, A. B., et al. 2019, ApJ, 885, 3

Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163

Sohn, S. T., Watkins, L. L., Fardal, M. A., et al. 2018, ApJ, 862, 52

Soubiran, C., Jasnowicz, G., Chemin, L., et al. 2013, A&A, 552, A64

Soubiran, C., Jasnowicz, G., Chemin, L., et al. 2018, A&A, 616, A7

Southworth, J., Gänsicke, B. T., Marsh, T. R., et al. 2006, MNRAS, 373, 687

Starkenburg, E., Martin, N., Youakim, K., et al. 2017, MNRAS, 471, 2587

Starkenburg, E., Youakim, K., Martin, N., et al. 2019, MNRAS, 490, 5757

Sun, G., Wang, Y., Liu, C., et al. 2023, Res. Astron. Astrophys., 23, 015013

Thomas, G. F., Famaey, B., Ibata, R., Lüghausen, F., & Kroupa, P. 2017, A&A, 603, A65

Thomas, G. F., McConnachie, A. W., Ibata, R. A., et al. 2018, MNRAS, 481, 5223

Thomas, G. F., Laporte, C. F. P., McConnachie, A. W., et al. 2019, MNRAS, 483, 3119

Thomas, G. F., Jensen, J., McConnachie, A., et al. 2020, ApJ, 902, 89

Titulaer, L. 2021, Bachelor's Thesis, University of Groningen, The Netherlands fse.studenttheses.ub.rug.nl/25625/

Tody, D. 1986, SPIE Conf. Ser., 627, 733

Tody, D. 1993, in Astronomical Society of the Pacific Conference Series, 52, Astronomical Data Analysis Software and Systems II, eds. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 173

Vasiliev, E. 2018, Astrophysics Source Code Library [record ascl:1805.008]

Vasiliev, E. 2019a, MNRAS, 482, 1525

Vasiliev, E. 2019b, MNRAS, 484, 2832

Vasiliev, E., & Baumgardt, H. 2021, MNRAS, 505, 5978

Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279

Vera-Ciro, C., & Helmi, A. 2013, ApJ, 773, L4

Vickers, J. J., Grebel, E. K., & Huxor, A. P. 2012, AJ, 143, 86

Vickers, J. J., Li, Z.-Y., Smith, M. C., & Shen, J. 2021, ApJ, 912, 32

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261

Viswanathan, A., Byström, A., Starkenburg, E., et al. 2024, A&A, submitted [arXiv:2408.17250]

Vivas, A. K., Zinn, R., & Gallart, C. 2005, AJ, 129, 189

Wang, J., Hammer, F., & Yang, Y. 2022, MNRAS, 510, 2242

Waskom, M. L. 2021, J. Open Source Softw., 6, 3021

Williams, A. A., & Evans, N. W. 2015, MNRAS, 454, 698

Williams, A. A., Belokurov, V., Casey, A. R., & Evans, N. W. 2017, MNRAS, 468, 2359

Xue, X. X., Rix, H. W., Zhao, G., et al. 2008, ApJ, 684, 1143

Xue, X.-X., Rix, H.-W., Yanny, B., et al. 2011, ApJ, 738, 79

Yanny, B., Newberg, H. J., Kent, S., et al. 2000, ApJ, 540, 825

Yanny, B., Rockosi, C., Newberg, H. J., et al. 2009, AJ, 137, 4377

York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579

Yu, F., Li, T. S., Speagle, J. S., et al. 2024, ApJ, 975, 81

Zhao, H. 1996, MNRAS, 278, 488

Zhou, Y., Li, X., Huang, Y., & Zhang, H. 2023, ApJ, 946, 73
## Appendix A: Reduction and preparation of observed and synthetic data

The systematic process of extracting the spectra from the two-dimensional data recorded by the MIT detector on the VLT/FORS2 instrument consists of four elements or components: (1) a basic correction for the artifacts and 'fingerprints' of the two-dimensional chip data, including bias level, flat field response, bad pixels, and slit illumination based on calibration data gathered at daytime (bias frames, flat fields) (2) a two-dimensional mapping (along the slit axis) of the dispersion axis to wavelengths with measured lines of known wavelength in arc-lamps at daytime and at zenith, (3) a flux calibration and instrument response correction from a response curve estimated from a measured spectrum in a set-up similar to that of the science observations of a flux standard star (usually within ±3 days of the observations of the science targets) with known continuum behaviour and fluxes, and (4) tracing the science spectrum along the slit in the chip data while adding together the measured pixel values at each mapped column within a defined aperture and subtracting the median sky value outside of the aperture along each CCD column. The calibration of the two-dimensional wavelength solution for the FORS2 observations with the 600B+22 grism included fitting a polynomial of fourth order with residuals typically of the order of 0.25 Å with a few outlier pixels at up to 1.00 Å, ensuring a wavelength calibration accuracy within typically 0.25 Å across the entire detector, allowing for an accurate representation of the spectral dispersion across the detector.

With respect to the behaviour of the spectral resolution across the detector for the 2D data taken with the slit with width of 1.0 arcsec (candidate BHB stars) and 0.4 arcsec (velocity standard stars), the version of the FORS2 pipeline within the ESOREFLEX that we employed, also provides an estimate of the average resolving power across the detector as measured from lines in the spectra of the arc-lamps across the dispersion axis for each wavelength pixel. From Fig. A.1 we can see that a line provides a sufficient description of the trends of the mean resolving power across the detector for each wavelength bin for both datasets and also various observation times together as we do not see any significant discrepancy in the trend between observation times (different markers in Fig. A.1). Utilising a $\chi^2$-minimisation using NUMPY.POLYFIT (Harris et al. 2020), we found the lines that best fit the observed trends in Fig. A.1 to be

$$R = 826 + 0.161 \AA^{-1} \times \lambda \qquad (A.1)$$

for slit width = 0.4 arcsec (left half of Fig. A.1) and

$$R = 293 + 0.119 \AA^{-1} \times \lambda \qquad (A.2)$$

as a function of wavelength $\lambda$ in Å for slit width = 1.0 arcsec (right half of Fig. A.1).

For the first observation of candBHB1 listed in Table 1 some calibration data—such as a normalised flat-field image in connection with the observations of a standard star to calibrate the flux of the spectrum of candBHB1—were missing and we were therefore not able to run the full data reduction and calibration of this observation with the pipeline. However, the first spectrum of candBHB1 was automatically reduced and calibrated by the ESO Science Data Quality Group in a set-up very similar to the one we adopted using the FORS2 pipeline 5.6.4, and we have adopted their results for this observation.

Regarding the synthetic spectra that we use to model the FORS2 spectra with RVSPECFIT, for the first instrumental set-up

connected to the spectra of the candidate BHB stars the synthetic stellar spectra of the PHOENIX library are re-binned to the dispersion axis of the FORS2 spectra between 3,900.94 Å and 5,989.18 Å (the wavelengths of the first and last wavelength bin after 3,900 Å and before 5,990 Å respectively, that we consider to be the reliable range of the spectra) with a linear dispersion = 1.32 Å and the spectral resolution was degraded from the initial $R = 100000$ to Equation A.2. All wavelengths are measured in air in this set-up.

Likewise, the library of synthetic stellar spectra that are used to fit the spectra of the velocity standard stars are prepared in the same way as the spectra of the candidate BHB stars with lower resolution except we have to adopt the synthetic spectra in a higher spectral resolution that behaves as Equation A.1.

Finally, the third grid of prepared synthetic stellar spectra that we use for the SDSS spectra have a linear dispersion of 1 Å across 3,500–9,500 Å and have a uniform resolving power of 2,000 across the dispersion axis (wavelengths are in vacuum).

For all three sets of synthetic spectra we perform an initial interpolation to produce a finer grid by using 0.2 dex as step size for both [Fe/H] and [α/Fe] while preserving the 0.5 dex sampling in $log_{10}(g)$ and non-uniform sampling in effective temperature. Similar interpolations have been used in the past in studies that also employ RVSPECFIT (e.g. Li et al. 2019; Ji et al. 2021; Cooper et al. 2023). After the first interpolation, the second interpolation is only performed during the fitting stage for each likelihood evaluation (cf. with comments in Li et al. 2019).

Table A.1 lists the definitions of the normally distributed priors of $T_{\text{eff}}$, $log_{10}(g)$, [Fe/H], and [α/Fe] that we adopted for the posterior calculation given within RVSPECFIT. These utilise estimates of $T_{\text{eff}}$, $log_{10}(g)$, [Fe/H], and [α/Fe] from Gaia DR3 the Gaia General Stellar Parametriser for Spectra catalogue (Gaia Collaboration 2023; Recio-Blanco et al. 2023) and sample statistics for the velocity standard stars for which no estimates of [Fe/H] and [α/Fe] was available in this catalogue.
<figure><img src="image_13.png" /><figcaption>**Fig. A.1.** VLT/FORS2/600B+22 grism resolving power trends across detector for gathered data of arc-lamps at different observation times for the two used slit configurations.</figcaption></figure>

<table><caption>**Table A.1.** Summary of normally distributed priors of RVSPECFit RVSPECFIT (Koposov et al. 2011; Koposov 2019) stellar parameters used for the inference of the LOS velocities from the VLT/FORS2/600B+22 grism spectra taken in the slit with width = 0.4 arcsec of the velocity standard stars with RVSPECFIT.</caption><thead><tr><th rowspan="2">Designation<br>(1)</th><th rowspan="2">T<sub>eff</sub><br>[K]<br>(2)</th><th rowspan="2">[K]<br>(3)</th><th colspan="2">log<sub>10</sub> g</th><th colspan="2">[Fe/H]</th><th colspan="2">[α/Fe]</th><th colspan="2">Source [Fe/H]</th><th colspan="2">Source [α/Fe]</th></tr><tr><th>(4)</th><th>(5)</th><th>(6)</th><th>(7)</th><th>(8)</th><th>(9)</th><th>(10)</th><th>(11)</th><th>(12)</th><th>(13)</th></tr></thead><tbody><tr><td>BD +24 1843</td><td>5570</td><td>33</td><td>4.17</td><td>0.11</td><td>-0.1</td><td>0.2</td><td>-0.01</td><td>0.0</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +31 1781</td><td>4532</td><td>36</td><td>4.37</td><td>0.05</td><td>-0.1</td><td>0.3</td><td>0.03</td><td>0.0</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +34 1955</td><td>5235</td><td>36</td><td>3.69</td><td>0.07</td><td>-0.4</td><td>0.3</td><td>0.04</td><td>0.0</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>WDS J07277+2420A</td><td>4201</td><td>150</td><td>4.24</td><td>0.20</td><td>-0.2</td><td>1.4</td><td>-0.2</td><td>0.1</td><td>med</td><td>max</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +26 1647</td><td>5894</td><td>510</td><td>4.89</td><td>2.55</td><td>-0.2</td><td>1.4</td><td>0.07</td><td>0.2</td><td>med</td><td>max</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +31 1684</td><td>4250</td><td>1500</td><td>1.89</td><td>2.84</td><td>-0.2</td><td>1.4</td><td>0.04</td><td>0.4</td><td>med</td><td>max</td><td>med</td><td>max</td></tr><tr><td>BD +30 1501</td><td>4540</td><td>33</td><td>2.11</td><td>0.08</td><td>-0.4</td><td>0.3</td><td>0.03</td><td>0.0</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>TYC 2461-988-1</td><td>6475</td><td>180</td><td>4.19</td><td>0.29</td><td>-0.2</td><td>1.4</td><td>0.2</td><td>0.4</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +36 1823</td><td>6041</td><td>45</td><td>3.98</td><td>0.05</td><td>-0.2</td><td>0.3</td><td>0.06</td><td>0.0</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr><tr><td>BD +35 1801</td><td>5972</td><td>66</td><td>4.09</td><td>0.11</td><td>-0.1</td><td>0.3</td><td>0.07</td><td>0.1</td><td>GSPS</td><td>GSPS</td><td>GSPS</td><td>GSPS</td></tr></tbody></table>

**Notes.** The mean of the Gaussian priors of the effective temperature (Col. (2)) are the estimates from the Gaia DR3 Gaia General Stellar Parametriser for Spectra (GSP-Spec) catalogue (Gaia Collaboration 2023; Recio-Blanco et al. 2023) that includes measurements from the Gaia Radial Velocity Spectrometer (Cropper et al. 2018). In Col. (3) the uncertainties of the GSP-Spec estimates of the effective temperature were scaled by the factor 3 to generate a broad Gaussian that mimics to some degree an uniform prior. In the same vein, we apply the same scaling for the other columns that give the standard deviations of the Gaussian priors (Col. (5), (7), and (9)). Similarly, the data in Col. (4) and (5) were derived from the GSP-Spec catalogue, using the calibration as recommended in Recio-Blanco et al. (2023) (equation 1 with polynomial coefficients in Recio-Blanco et al. 2023) to correct the bias for their estimates. The means and standard deviations of the broad Gaussian priors of [Fe/H] (Col. (6) and (7)) and [α/Fe] (Col. (8) and (9)) are based on calibrated values from the GSP-Spec catalogue if available following the recommendations in Recio-Blanco et al. (2023), that is, in case of [M/H] using equation 2 with polynomial coefficients given in second row of table 3 in Recio-Blanco et al. (2023) to calculate [Fe/H] = [Fe/M] + [M/H] and with respect to [α/Fe] using equation 3 with coefficients of the fourth-order polynomial in table 4 in Recio-Blanco et al. (2023) (second [α/Fe] of table 4). For some of the velocity standard stars, no GSP-Spec estimates of [M/H] and/or [α/Fe] were available. For these stars we took the median of all the other means of [Fe/H] and [α/Fe] as the mean of the Gaussian prior and the maximum of the other standard deviation entries regarding the standard deviation of the Gaussian prior. In the last four columns GSPS is used here to refer to estimates from the GSP-Spec catalogue, 'med' indicates that for the entry of the mean of the Gaussian prior of [Fe/H] or [α/Fe] we used the median approach as explained above, and 'max' labels an entry where we adopted the maximum standard deviation of all Gaussian priors of [Fe/H] or [α/Fe] as also noted above.

## Appendix B: Pan-STARRS1 (g-r)₀-effective temperature relation of (candidate) blue horizontal branch stars

Gaussian priors of effective temperature were selected for the posterior calculation given the spectra of the candidate BHB stars in the sample. For the estimation of effective temperature of each candidate BHB star in the sample, a relation between (g-r)ₚ₅₀ and effective temperature for the BHB stars in the Barbosa et al. (2022) sample (probability of being BHB star of 100 per cent) was derived. We adopt PS1 DR2 mean forced photometry (Chambers et al. 2016). Estimates of effective temperature from Sloan Extension of Galactic Understanding and Explo-

ration (SEGUE) Stellar Parameter Pipeline (Lee et al. 2008a,b) based on SEGUE data (Yanny et al. 2009) of BHB stars were chosen because data on these estimates are already included in the Barbosa et al. (2022) catalogue. An orthogonal distance regression was used to find the best-fitting linear relation between (g-r)ₚ₅₀ and effective temperature for the BHB stars in the set released by Barbosa et al. (2022) and also allow to take into account uncertainties for both quantities in fitting. For the fitting we used the SCIPY.ODR package (Virtanen et al. 2020). The results of the correlational analysis are summarised in Fig. B.1. In Fig. B.1 there is a clear trend of decreasing effective temperature with increasing (g-r)ₚ₅₀ that can be described by
<figure><img src="image_14.png" /><figcaption>**Fig. B.1.** Correlation of $(g-r)_{PS,0}$ and effective temperature among BHB stars in the Barbosa et al. (2022) sample (probability of being BHB star of 100 per cent)</figcaption></figure>

$T_{\text{eff}} = -(5349 \pm 54) \text{ K} \times (g-r)_{\text{PS,0}} + (7504 \pm 8) \text{ K}$ and residual variance of the order of 5 K (line in Fig. B.1).

## Appendix C: Comparison of line-of-sight velocities of velocity standard stars

Our comparison of the maximum a posteriori estimates of the heliocentric LOS velocities of the velocity standard stars can be seen in Table 2 from fitting the optical low-resolution spectra taken with the FORS2 instrument using the 600B+22 grism in a 0.4 arcsec wide slit with RVSPECFIT and reference values by Soubiran et al. (2018). Further analysis suggests a slight drop of the difference with higher median signal-to-noise ratio of the spectra between 3,900Å and 5,990Å (considered the reliable range of the spectra) in the upper panel of Fig. C.1 and no evidence of any relational association between the LOS velocity difference and mid-observation time (lower panel of Fig. C.1). Overall, these results indicate that it is sufficient to investigate this comparison as a function of air mass under which the stars were observed as we do in Sect. 2.6.

## Appendix D: Sample of candidate M-type giant stars

A sample of (candidate) M-type giants identified by using 822,752 low-resolution spectra from LAMOST DR9 (Cui et al. 2012) that cover the spur of the Sagittarius stream (Li et al. 2023a). Li et al. (2023a) present a catalogue of their (candidate) M-type giants. (Li et al. 2023a) reported 183 candidate members of the Sagittarius stream.

In order to compile their catalogue of (variable) RR Lyrae stars in the Sagittarius stream, Hernitschek et al. (2017) consider all RRab stars identified in Sesar et al. (2017b) to be part of the Sagittarius stream if they are within nine degrees of the Sagittarius stream orbital plane. We have adopted the same criteria for the selection of candidate Sagittarius M-type giant stars in the Li et al. (2023a) catalogue as we for the selection of RR Lyrae in Fig. 5. The positions of all stars in the Li et al. (2023a) catalogue on the Sagittarius stream orbital plane, as defined in Vasiliev et al. (2021), were calculated using the GALA (Price-Whelan 2017) software. Quantitatively, the Hernitschek et al. (2017) selection criterion can be expressed as $|B_{Sgr}| < 9^\circ$ where

<figure><img src="image_15.png" /><figcaption>**Fig. C.1.** Results obtained from the analysis of the FORS2 spectra of the velocity standard stars compared to the reference LOS velocity values by Soubiran et al. (2018). The times of mid-observations of the spectra are shown in the abscissa of the second panel.</figcaption></figure>

$B_{Sgr}$ is the latitude of the Vasiliev et al. (2021) definition of the Sagittarius stream coordinate system.

Several M-giants are observed in the region of the sky of the spur region and consistent with its distance range. However, the catalogue of Li et al. (2023a) provides no uncertainties of their photometric distances. Their distances are estimated from a empirical Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006) colour $(J-K)_0$-absolute magnitude $M_{J0}$ relation (equa- tion 3 in Li et al. 2023a) fitted with a subset of the total sample of (candidate) M-type giants closer than 4 kpc and Gaia early DR3 inverted relative parallax uncertainties (Gaia Collaboration 2021) larger than 5 in the Bailer-Jones et al. (2021) catalogue. However, in their figure 5, Li et al. (2023a) examine the extent to which the calculated $M_{2\text{MASS J0}}$ can be trusted and find a spread of the relation of 0.64 dex. This residual scatter of 0.64 will be presumably the largest component in the uncertainty of the distance modulus in 2MASS *J* that they use to compute the pho- tometric heliocentric distance of the (candidate) M-type giants in the catalogue.

An uncertainty of 0.64 dex would translate to a photometric $\Delta D_\odot$ of at least 30 kpc for stars at a distance of $\ge 70$ kpc as we have in our sample. Therefore, we choose not to use these stars for our purposes.