#  ${S}^{5}$: Probing the Milky Way and Magellanic Clouds potentials with the
6-D map of the Orphan-Chenab stream

Sergey E. Koposov,1,2,3, Denis Erkal4, Ting S. Li5, Gary S. Da Costa6,7, Lara
R. Cullinane8, Alexander P. Ji9,10, Kyler Kuehn11,12, Geraint F. Lewis13,
Andrew B. Pace14, Nora Shipp15, Daniel B. Zucker16,17, Joss Bland-
Hawthorn13,7, Sophia Lilleengen4, and Sarah L. Martell18,7(${S}^{5}$
Collaboration)
1 Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK
2 Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge
CB3 0HA, UK
3 Kavli Institute for Cosmology, University of Cambridge, Madingley Road,
Cambridge CB3 0HA, UK
4 Department of Physics, University of Surrey, Guildford GU2 7XH, UK
5 Department of Astronomy and Astrophysics, University of Toronto, 50 St.
George Street, Toronto ON, M5S 3H4, Canada
6 Research School of Astronomy and Astrophysics, Australian National
University, Canberra, ACT 2611, Australia
7 Centre of Excellence for All-Sky Astrophysics in Three Dimensions (ASTRO
3D), Australia
8 Department of Physics and Astronomy, Johns Hopkins University, 3400 N.
Charles St, Baltimore, MD 21218, USA
9 Department of Astronomy & Astrophysics, University of Chicago, 5640 S Ellis
Avenue, Chicago, IL 60637, USA
10 Kavli Institute for Cosmological Physics, University of Chicago, Chicago,
IL 60637, USA
11 Lowell Observatory, 1400 W Mars Hill Rd, Flagstaff, AZ 86001, USA
12 Australian Astronomical Optics, Faculty of Science and Engineering,
Macquarie University, Macquarie Park, NSW 2113, Australia
13 Sydney Institute for Astronomy, School of Physics, A28, The University of
Sydney, NSW 2006, Australia
14 McWilliams Center for Cosmology, Carnegie Mellon University, 5000 Forbes
Ave, Pittsburgh, PA 15213, USA
15 MIT Kavli Institute for Astrophysics and Space Research, 77 Massachusetts
Ave., Cambridge, MA 02139, USA
16 School of Mathematical and Physical Sciences, Macquarie University, Sydney,
NSW 2109, Australia
17 Macquarie University Research Centre for Astronomy, Astrophysics &
Astrophotonics, Sydney, NSW 2109, Australia
18 School of Physics, UNSW, Sydney, NSW 2052, Australia
E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present a 6-D map of the Orphan-Chenab (OC) stream by combining the data
from Southern Stellar Stream Spectroscopic Survey (${S}^{5}$) and Gaia. We
reconstruct the proper motion, radial velocity, distance, on-sky track and
stellar density along the stream with spline models. The stream has a total
luminosity of $M_{V}=-8.2$ and metallicity of $\mathrm{[Fe/H]}=-1.9$, similar
to classical Milky Way (MW) satellites like Draco. The stream shows drastic
changes in its physical width varying from 200 pc to 1 kpc, but a constant
line of sight velocity dispersion of 5 $\mathrm{\,km\,s^{-1}}$. Despite the
large apparent variation in the stellar number density along the stream, the
flow rate of stars along the stream is remarkably constant. We model the 6-D
stream track by a Lagrange-point stripping method with a flexible MW potential
in the presence of a moving extended Large Magellanic Cloud (LMC). This allows
us to constrain the mass profile of the MW within the distance range 15.6 < r
< 55.5 kpc, with the best measured enclosed mass of $(2.85\pm 0.1)\times
10^{11}\,\mathrm{\,M_{\odot}}$ within 32.4 kpc. Our stream measurements are
highly sensitive to the LMC mass profile with the most precise measurement of
its enclosed mass made at 32.8 kpc, $(7.02\pm 0.9)\times 10^{10}\,{\rm
M}_{\odot}$. We also detect that the LMC dark matter halo extends to at least
53 kpc. The fitting of the OC stream allows us to constrain the past LMC
trajectory and the degree of dynamical friction it experienced. We demonstrate
that the stars in the OC stream show large energy and angular momentum spreads
caused by LMC perturbation.

###### keywords:

Key words: Galaxy: evolution – Galaxy: halo – Galaxy: kinematics and dynamics
– Galaxy: structure – Magellanic Clouds

††pubyear: 2022††pagerange: ${S}^{5}$: Probing the Milky Way and Magellanic
Clouds potentials with the 6-D map of the Orphan-Chenab stream–28

##  1 Introduction

Stellar streams are one of the most impressive illustrations of the
hierarchical nature of galaxy formation predicted by the $\Lambda CDM$
paradigm (White & Frenk, 1991). While Milky Way-like galaxies are expected to
have accreted a large number of stellar systems over their assembly history
(Bullock & Johnston, 2005), only in the last 20 years or so have we started to
obtain direct evidence of a large number of disrupting low mass systems in the
Milky Way stellar halo (e.g. Ibata et al., 1994; Odenkirchen et al., 2003;
Belokurov et al., 2006; Grillmair & Dionatos, 2006).

Since the discovery of the first streams in the Milky Way halo, it was clear
that streams can provide powerful constraints on the Galactic potential
(Johnston et al., 1999; Binney, 2008; Koposov et al., 2010) as streams
approximately trace the orbit of their progenitor. While initial models of
streams were often quite simplistic and relied on the (incorrect) assumption
that streams are orbits (Sanders & Binney, 2013), new, more sophisticated ways
of modelling stellar streams have since been developed (Bovy, 2014; Gibbons et
al., 2014; Fardal et al., 2015). With the number of known streams growing over
time due to increased depth and coverage of Galactic surveys (Koposov et al.,
2014; Bernard et al., 2014), this has opened a window for methods that fit
multiple stellar streams together with the aim of breaking degeneracies
between the parameters of the Galactic potential (Bovy et al., 2016; Bonaca &
Hogg, 2018). However, as yet more streams were discovered with the DES survey
(Shipp et al., 2018) and Gaia (Malhan et al., 2018), the idealistic picture of
streams being a clean and easy tracer that can be used to constrain the local
Galactic potential has become more muddied.

As streams were mapped in greater detail using photometric data, it has become
clear that streams do not look like clean Gaussian trails of stars. Instead
they show stream gaps (Carlberg & Grillmair, 2013; Erkal et al., 2017),
broadenings (Price-Whelan et al., 2016; Malhan et al., 2018), bifurcations
(Belokurov et al., 2006), and "wiggles" (Li et al., 2021). The exact
attribution of these structures to specific physical effects is still a matter
of debate but is likely some combination of perturbations by other MW
satellites (de Boer et al., 2020; Li et al., 2021; Dillamore et al., 2022), DM
subhalos (Ibata et al., 2002; Johnston et al., 2002; Yoon et al., 2011;
Carlberg, 2020), giant molecular clouds (Amorisco et al., 2016), and
interaction with the Milky Way bar (Hattori et al., 2016). Together this makes
the analysis of streams and their use as a confident tracer of the MW
potential more difficult.

The final nail in the coffin of streams being straightforward tracers of the
MW potential was the study of the so-called Orphan stream (Grillmair &
Dionatos, 2006; Belokurov et al., 2007) with Gaia DR2 data. While initially
the stream was traced over approximately 60 degrees on the sky, the RR Lyrae
from Gaia DR2 allowed Koposov et al. (2019) (hereafter K19) to trace it for
$\sim 200$ degrees. It turned out that the stream extension to the South was
actually previously detected and was called "Chenab" by Shipp et al. (2018).
The reason for the mis-association of the stream as a new stream, as opposed
to the continuation of Orphan, was that their orbital planes were seen as
quite distinct. However, the cause for this was found to be that the stream is
being actively perturbed by the Large Magellanic Cloud (Erkal et al. 2019,
hereafter E19), causing the stream to twist and the stars in the stream to
move sideways. While complicating the analysis of the Milky Way potential,
this nevertheless allowed them to put strong constraints on the mass of the
Large Magellanic Cloud (LMC), revising its mass to about a fifth or tenth of
the Milky Way’s. This in turn led to the realisation of a very strong
perturbation to the Milky Way potential by the Clouds that manifests itself in
multiple forms: direct gravitational pull by the Clouds and second-order
effects such as the LMC DM wake (Garavito-Camargo et al., 2019; Belokurov et
al., 2019; Conroy et al., 2021), the motion of the MW centre around the LMC+MW
centre of mass (Gómez et al., 2015; Erkal et al., 2021; Petersen & Peñarrubia,
2021), and the deformation of the MW dark matter distribution (Lilleengen et
al., 2023). Since then, the effects of the LMC, expressed as stars moving at
an angle with respect to the stream itself, have been seen in other streams
(Shipp et al., 2019).

The effects of the LMC have made stream modelling a much more difficult task,
since it requires introducing a time-dependent LMC potential, and accounting
for the motion of the MW centre. However, this modelling should now allow us
to start probing not only the MW dark matter distribution, but also the LMC
dark matter distribution, and how these deform throughout the interaction.
Such sophisticated models require extensive, ideally fully 6-D datasets. Erkal
et al. (2019) developed the first such model to fit the 5-D data of the OC
stream from Koposov et al. (2019). Vasiliev et al. (2021) have also
demonstrated this with a model of the Sagittarius stream, where model took
into account the MW and LMC interaction to constrain the parameters of both MW
and LMC.

In this paper, we focus on the Orphan-Chenab (or OC) stream as it is probably
the best stream for probing the MW potential and the LMC. This stream is very
long, probes the Galaxy from $\sim 15$ kpc to $\sim 80$ kpc, and is affected
by the LMC. It also does not have the complexity of the Sgr stream with its
unexplained bifurcation. Because of this, the OC stream was followed up by the
Southern Stellar Stream Spectroscopic Survey (${S}^{5}$) over the last 5 years
in order to make a detailed phase space map of the stream (Li et al., 2019).
In this paper, we combine the results of the dedicated ${S}^{5}$ follow-up of
the OC stream with other public spectroscopic datasets and recent Gaia data to
model the 6-D phase space of the OC stream. We then use this 6-D track to
infer the potential of the Milky Way and LMC through a detailed stream model
fitted to the 6-D stream track.

In the paper we also discuss possible origins of the OC stream and location of
the stream progenitor, as despite many claims in the literature of possible
association of the stream with known objects (Fellhauer et al., 2007;
Grillmair et al., 2015; Koposov et al., 2019), no stream progenitor has been
identified yet.

The paper is structured as follows. In Section 2 we describe the data we rely
on in this paper. Section 3 demonstrates how we model the 6-D phase space
track of the stream. In Section 4 we analyse and interpret observational
features extracted in Section 3. Section 5 describes the stream model we
employ and in Section 6 we discuss the inference from the model. We conclude
in Section 7.

##  2 Data

###  2.1 Observations and Target Selection for ${S}^{5}$

Most of the data presented in this paper are collected as part of the
observations by ${S}^{5}$, which couples the Two-degree Field (2dF) fiber
positioner (Lewis et al., 2002) with the dual-arm AAOmega spectrograph (Sharp
et al., 2006) on the 3.9-m Anglo-Australian Telescope (AAT) to pursue a
complete census of known streams in the Southern Hemisphere. For more details
about observation, instrument setup and target selection, we refer readers to
Li et al. (2019). The OC stream is one of the 20 streams that have been mapped
by ${S}^{5}$.

This paper includes ${S}^{5}$ observations between August 2018 and April 2021.
A total of 41 fields dedicated to the OC stream were observed at the AAT, 34
of which were observed in 2018 and 2019, and are listed in Table 2 of (Li et
al., 2019); the remaining 7 were added in 2020 and 2021. All of the fields
have an exposure time between 3000 and 7200 seconds. Since the OC stream
covers a large area on the sky and some other streams intersect it, we also
include these stream fields in addition to the dedicated OC fields.

As described in Li et al. (2019), we identify probable OC member stars by
applying selections in proper motion space and in colour-magnitude diagram
(CMD) space. In particular, we select the targets based on the proper motions
and distances of the RR Lyrae stars in OC from K19. We compute the difference
in proper motions between the targets and the average proper motion from RR
Lyrae stars at a given stream longitude, or $|\Delta\mathrm{PM}|$. For the
majority of fields we used $|\Delta\mathrm{PM}|<2$ $\mathrm{\,mas\,yr^{-1}}$,
while close to the Galactic plane (at $|\phi_{1}|\lesssim 15$ deg) we used
$|\Delta\mathrm{PM}|<1$ $\mathrm{\,mas\,yr^{-1}}$. For most of the fields DES
photometry was used for CMD selection; otherwise, Gaia photometry was used.

###  2.2 Data Reduction for ${S}^{5}$

In this paper, we rely on the ${S}^{5}$ data from the internal release DR3.1
which differs from the one described in Li et al. (2019) and released in Li &
S5 Collaboration (2021) in several aspects. We have significantly larger
coverage of the sky compared to DR1, as well as improved modelling. For more
details of the model fitting updates, we refer readers to Li et al. (2022). In
particular, we are now modelling blue and red spectra of all the exposures
simultaneously for each object in order to improve stellar parameter
estimates. The DR3.1 data also provides the spectro-photometric distance
estimates, which are obtained by simultaneous fitting of the spectroscopy,
photometry and Gaia parallaxes by interpolated PHOENIX (Husser et al., 2013)
spectral templates and MIST isochrones (Koposov et al. 2023 in prep).

###  2.3 Other spectroscopic surveys

In this paper, we use radial velocity and iron abundance measurements from
${S}^{5}$ and from several large spectroscopic surveys: APOGEE DR17 (Nidever
et al., 2015; Ahumada et al., 2020; Jönsson et al., 2020), LAMOST DR7 (Cui et
al., 2012; Zhao et al., 2012), and SDSS DR14 (Abolfathi et al., 2018)
measurements by the SSPP spectroscopic pipeline (Lee et al., 2008a, b;
Smolinski et al., 2011). We use the vhelio_avg, verr and fe_h columns from
allstar table in APOGEE, and rv, rv_err, and feh columns from LAMOST LRS
table111 We apply a 5 $\mathrm{\,km\,s^{-1}}$ offset to LAMOST radial
velocities as this is the typical velocity offset we notice with respect to
APOGEE velocities.. We use the fehadop, elodiervfinal, and elodiervfinalerr
columns from SDSS.

###  2.4 Gaia

Throughout the paper, unless specified otherwise, we rely on Gaia  EDR3 data
(Gaia Collaboration et al., 2021), such as parallaxes, proper motions and
BP/RP/G magnitudes.

##  3 6-D track extraction

In this section, we describe the modelling of the OC stream stars’ positions,
velocities, distances and colour-magnitude diagram locations in order to
extract the 6-D phase space track of the stream.

Figure 1: Identification of spectroscopic OC stream members in the ${S}^{5}$,
SDSS, LAMOST and APOGEE datasets. The 4 panels show the distribution of
various properties of stars in our sample (position on the sky $\phi_{2}$, GSR
radial velocity, and proper motion in $\phi_{1}$ and $\phi_{2}$) versus the
angle along the OC stream ($\phi_{1}$). The stars plotted in each panel have
been selected using all other panels: i.e. the plot of $\phi_{2}$ vs
$\phi_{1}$ relies on proper motion and radial velocity selection but not on
position. The grey lines delineate the selection boundaries for likely OC
stream members. Black points show stars from ${S}^{5}$, red from SDSS, blue
from LAMOST and green from APOGEE. The stream is excellently traced throughout
200 degrees on the sky (c.f. Figure 2 from K19, ).

###  3.1 Conventions

Throughout the paper, we rely on the standard coordinate system aligned with
the Orphan stream ($\phi_{1}$, $\phi_{2}$) as defined in K19. As the
$(\phi_{1},\phi_{2})=(0,0)$ point is the location where the stream crosses the
MW plane and $\phi_{1}>0$ part of the stream is in the North Galactic Cap, we
will refer the $\phi_{1}>0$ parts of the stream as Northern, while
$\phi_{0}<0$ as Southern. We also use $\mu_{\phi,1}$, $\mu_{\phi,2}$ which
refer to the proper motions in the $\phi_{1}$, $\phi_{2}$ heliocentric frame
(not corrected for the solar motion). The proper motions in right ascension
and $\phi_{1}$ (unless specified otherwise) always include the cosine term.

In the paper, we occasionally use measurements in the Galactocentric frame.
This frame is based on the astropy v4.0 frame (Reid & Brunthaler, 2004;
Gravity Collaboration et al., 2018; Drimmel & Poggio, 2018; Bennett & Bovy,
2019).

###  3.2 Kinematic sample

We proceed to map the OC stream in steps with the goal of fully extracting its
track in 6-D phase space. We start with the spectroscopic data as these
provide us with the cleanest separation of stream members from foreground
stars allowing us to bootstrap our OC member identification.

We first assemble the list of potential OC members from the ${S}^{5}$, APOGEE,
SDSS, and LAMOST surveys by selecting the stars within 10 degrees of the great
circle established by K19. On top of the spatial selection, we apply a Gaia
parallax selection, $\omega<3\sigma_{\omega}+0.1$ mas, which should remove
nearby (closer than $\sim$ 10 kpc) stars. The small number of duplicates
observed in multiple surveys are removed. When deciding which duplicate to
remove we prioritise the surveys such that highest priority is given to APOGEE
followed by ${S}^{5}$, SDSS and LAMOST. We require that the spectroscopic
metallicity is $\mathrm{\,[Fe/H]}<-1.5$. For the ${S}^{5}$ stars we also
require that the 84-th percentile of spectro-photometric distance is larger
than 10 kpc (this guarantees that the star is confidently not a nearby star)
to further reduce the contamination from nearby stars. Figure 1 shows the
spatial and kinematic distribution of the sample as a function of the angle
along the stream. On the plot we also show the linear spline-based regions
(grey lines) that we use to make the preliminary selection of possible OC
stream members. We call these fiducial splines, and they are broadly based on
the stream tracks from K19. The splines are provided in the supplementary
materials of the paper and the width of the selection region is 3 degrees for
the spatial selection, 30 $\mathrm{\,km\,s^{-1}}$ for the radial velocity
selection, and 0.5 $\mathrm{\,mas\,yr^{-1}}$ in proper motion. The total
number of possible stream members from different surveys that fall within the
spline boundaries in the $-100<\phi_{1}<150$ range is 213, 110, 25 and 12 for
${S}^{5}$, SDSS, LAMOST and APOGEE respectively.

For Figure 1, we use a convention similar to the one adopted in Figure 2 of
K19, where for each panel of the plot we use the selection from all the other
panels. In other words, the plot of radial velocity versus angle only shows
stars selected based on positions and proper motions (but not radial
velocity), while the plot of the $\mu_{\phi,1}$ proper motion versus angle
includes stars selected based on $RV$ and position on the sky and
$\mu_{\phi,2}$ (but not $\mu_{\phi,1}$). This type of plot allows us to
cleanly see the completeness and contamination of our selections, as well as
better assess how prominent the stream signal is. The coloured points on the
plot show measurements from different surveys, with black points representing
the ${S}^{5}$ data. The filled light grey regions in the top panel show areas
targeted by the ${S}^{5}$ survey.

The OC stream track on Figure 1 is clearly visible throughout almost 200
degrees on the sky from $\phi_{1}\sim-80$ deg to $\phi_{1}\sim 120$ deg on all
panels, matching the behaviour tracked by RR Lyrae in K19, but with
significantly more stars, higher proper motion precision, and the measurement
of radial velocities. We remark that the ${S}^{5}$ data alone maps the stream
over almost $\sim 110$ degrees without any significant gaps, while for the
Northern part of the stream ($\phi_{1}\gtrsim 80$ deg), the ${S}^{5}$ data is
complemented by SDSS, LAMOST, and somewhat by APOGEE.

In the next several sections, we will now rely on either the spectroscopic
sample of likely stream member stars constructed with the splines boundaries
shown on Figure 1, or the same splines selections applied to Gaia-only data
without spectroscopy. The spectroscopic sample of possible OC members based on
fiducial splines shown on Figure 1 is provided in Table 3.

###  3.3 RR Lyrae Distance model

In this section, we construct a model of the distance along the stream. We use
the catalogue of RR Lyrae from Gaia DR2 (Holl et al., 2018; Clementini et al.,
2019) complemented by the RR Lyrae from Stringer et al. (2021) and Sesar et
al. (2017) cross-matched with Gaia EDR3. We compute distances from Gaia
magnitudes using the equation from Table 4 of Muraveva et al. (2018) and the
extinction coefficient of $A_{G}/E(B-V)=2.27$ from Iorio & Belokurov (2021).
We use the phot_g_mean_mag average magnitudes from Gaia EDR3, shifted by 0.03
mag due to the typical offset between phot_g_mean_mag and the luminosity
averaged mean magnitude used by Muraveva et al. (2018).

Figure 2: Modelling of the distances to the OC stream stars. Top panel: The
distance moduli to the RR Lyrae selected by proper motion and position on the
sky. Grey lines show the selection boundaries for stars considered in the fit.
Middle panel: Grey points show the distance moduli of RR Lyrae selected by
proper motion and position on the sky. Black points with error-bars show the
distance modulus spline knots. Red lines show sampled models from the
posterior. Bottom panel: Base 10 logarithm of intrinsic distance modulus
spread at the location of spline knots inferred from the model. We observe an
increase of the distance modulus spread to $\sim$ 0.2 mag in the North from
about 0.05 mag in the South.

For the list of assembled RR Lyrae, we compute the positions and proper
motions along $\phi_{1}$ and $\phi_{2}$, selecting only the ones where
$\mu_{\phi,1}$ and $\mu_{\phi,2}$ are within 0.25 mas yr-1 of the fiducial
proper motion tracks shown Figure 1 and within 2 degrees of the stream track
on the sky. We then describe the stream distance modulus versus angle along
the stream using the spline-based mixture model:

$\displaystyle
P(m-M|\phi_{1})=(1-f(\phi_{1}))U(m-M|{\mathcal{D}}_{l}(\phi_{1}),{\mathcal{D}}_{h}(\phi_{1}))+$
$\displaystyle
f(\phi_{1})N^{*}(m-M|{\mathcal{D}}(\phi_{1}),\mathcal{S_{\mathcal{D}}}(\phi_{1}))$
(1)

Here, we model the distance modulus distribution as a mixture model of the
stream stars and contamination, where the stream stars are modelled as a
truncated Gaussian ($N^{*}$) distributed around the stream track
${\mathcal{D}}(\phi_{1})$ and with the $\phi_{1}$-dependent width
$\mathcal{S_{\mathcal{D}}}(\phi_{1})$, both of which are represented by cubic
splines. We model the data with distance moduli between
${\mathcal{D}}_{l}(\phi_{1})$, and ${\mathcal{D}}_{h}(\phi_{1})$ which is a 3
mag wide interval roughly centred on the distance track from K19. The
background is modelled as a uniform distribution, and the mixing fraction
between the background and the stream $f(\phi_{1})$ is modelled by a
spline222Here, and for models defined later, since the cubic spline cannot be
made bounded between zero and 1, for the mixing fraction spline we use the
transformation from the real line into [0,1] range using either
$F(x)=\frac{1}{2}+\frac{1}{\pi}\arctan x$ or logit function (depending on the
model) and add an appropriate Jacobian to the posterior (see e.g. K19, ).. For
this and all spline-based models described later we typically use different
number of spline knots for different parameters. The number of knots is chosen
by hand, motivated by constraining power of data and quality of the fit. The
model is implemented using the Stan probabilistic programming language
(Carpenter et al., 2017) and is sampled using the CMDStan and CMDstanpy
packages. The splines were implemented using the stan-splines package
(Koposov, 2022a)333https://github.com/segasai/stan-splines. The Stan code of
this and other models used in the paper is provided in supplementary materials
(see Data availability section). 444Here and for all the Stan models described
in the next sections, we always run the posterior sampling for at least 10000
iterations using a number of parallel chains (between 6 and 36), and ensure
that the convergence statistic $\hat{R}$ is below 1.01..

The results of the model are shown in Figure 2 and are also given in Tables 4
and 5 in the Appendix. We note that our model does not have spline knots at
$\phi_{1}=0$ deg, which is the point where the stream crosses the disc plane.
Our model is not constrained there as we do not have likely RR Lyrae members
between $-10\,{\rm deg}<\phi_{1}<20$ deg. This lack of constraints can be seen
in the spread of red curves at $\phi_{1}=0$ deg in the middle panel of Figure
2. We remark that there is a noticeable increase in the inferred spread of the
distance modulus in the North at $\phi_{1}>60$ deg which corresponds to a
$\sim$ 5% (or $\sim 4\,$kpc) distance spread that was previously noted by
Sesar et al. (2013). The spread can also be visually seen in the top panel of
the Figure (especially when compared to a small spread in the South at
$\phi_{1}\sim-50$ deg).

###  3.4 Colour magnitude diagram

We can use the sample of stars identified in Figure 1 to make a colour-
magnitude diagram of likely spectroscopic stream members. We select all stars
falling within the spline boundaries shown in the Figure, and use those to
make a Gaia-based colour-magnitude diagram (CMD). We correct the BP, RP and G
magnitudes by the distance modulus model derived in the previous section as
well as for extinction. Figure 3 shows the resulting CMD with points
representing stars selected to be within the position, proper motion, radial
velocity bounds on Figure 1, and the colour of each point indicating the
associated survey. Grey points show the stars that have positions and proper
motions within the selection boundaries but are not within the radial velocity
boundary of 30 km s-1. The spectroscopic members clearly form a well-defined
red giant branch and a horizontal branch with an RR Lyrae gap, with only a
handful of stars lying far from these branches. The grey points on other hand
are very diffusely distributed on the plot, indicating that our sample is very
pure. Although we can see stream members on the giant and subgiant branches,
the main sequence turn-off and the main sequence are not probed due to the
typical distance to OC stream members of more than 20 kpc and spectroscopic
magnitude limits.

To enable an efficient OC candidate member selection, in the next sections, we
use Figure 3 to define a polygon (shown in grey) that encompasses the OC
stellar populations in CMD space. We remark that we decided to avoid using the
model isochrones to perform CMD selection as we have struggled to find one
amongst several isochrone libraries for the Gaia photometry that well
describes both the RGB and the horizontal branch, and has a metallicity
matching our spectroscopic measurements.

Figure 3: The colour - absolute magnitude diagram of likely OC stream members.
The stars have been shifted by the best-fit distance model from Section 3.3.
Stars that fall within the proper motion, position and radial velocity
selection from Figure 1 are shown by points of different colours indicating
the origin of the spectroscopy. The greyed points are stars outside the
selection region in radial velocity space. The grey contour shows the adopted
colour magnitude selection of likely stream stars.

###  3.5 Radial velocity modelling

Figure 4: The measurement of radial velocities of the OC stream. Top panel:
The Galactic Standard of rest velocities of the stream as a function of angle
along the stream. The stars were selected by proper motions, position on the
sky, CMD and spectroscopic metallicity. Middle panel: The measurement of
velocities at spline knots is shown by black points with error-bars. The red
curves show the samples from the posterior. The grey points show the fitted
sample of stars. Bottom panel: The velocity dispersion of the OC stream as
measured by the spline model. The error-bars show the measurements at the
spline knots.  Figure 5: OC stream stars proper motions and our models. Top
panel: Black points show the measured (not-reflex corrected) proper motions
along $\phi_{1}$ vs angle along the stream. The stars were selected by CMD,
radial velocities, and stream track on the sky. Grey lines show the range of
proper motions that were modelled. The red curves show the samples from the
posterior. Second panel: The measurements of intrinsic (i.e. after taking into
account measurement errors) proper motion dispersion in $\mu_{\phi,1}$ from
the model. Third panel: The intrinsic dispersion of proper motions converted
to physical velocity. Bottom three panels: Same as first three panel but for
$\mu_{\phi,2}$.

In this section, we extract the radial velocity track of the OC stream. We use
stars with spectroscopic data identified in Section 3.2. We apply the proper
motion selection, colour-magnitude diagram selection and stream track on the
sky selection. We then construct the model for the radial velocity variation
along the stream. We model the radial velocities in the Galactic frame as a
mixture between the stream, represented by a Gaussian with variable mean and
velocity dispersion, and contamination (from the stellar halo) represented by
a Gaussian. The likelihood function for the radial velocity conditional on
$\phi_{1}$ is given below:

$\displaystyle
P(v|\phi_{1})=(1-f(\phi_{1}))N(v|v_{bg}+d_{bg}\phi_{1},{\mathcal{S}}_{bg})+$
$\displaystyle
f(\phi_{1})N(v|{\mathcal{V}}(\phi_{1}),{\mathcal{S}}_{v}(\phi_{1}))$ (2)

Here, $f(\phi_{1})$ is the mixing fraction, ${\mathcal{V}}(\phi_{1})$ is the
velocity of the stream, and ${\mathcal{S}}_{v}(\phi_{1})$ is the velocity
dispersion of the stream, all of which are represented as natural cubic
splines with the values at the knots being model parameters. $v_{bg}$ is the
centre of the Gaussian for the background component, $d_{bg}$ is the gradient
of the mean velocity of the background with respect to $\phi_{1}$, and
${\mathcal{S}}_{bg}$ is the velocity dispersion for the background. This model
(as other models in the paper) is implemented in the Stan programming language
and the code is provided in the supplementary materials. When computing the
likelihoods of individual stars, we account for their radial velocity
uncertainties by adding them in quadrature to the intrinsic velocity
dispersions. The priors on model parameters are mostly non-informative, where
the prior for mixing fraction at spline knots is $U(0,1)$, $N(0,300)$ for the
stream velocities, and $N(1.6,1.6)$ for the natural logarithm of velocity
dispersion.

The results of the modelling are shown in Figure 4. The figure is structured
similarly to Figure 2. The top panel shows the radial velocities in the
Galactic frame as a function of the angle along the stream for the modelled
subset of stars. The middle panel shows the knots of the RV spline with their
corresponding error-bars, as well as samples from the posterior for the
${\mathcal{V}}(\phi_{1})$ spline curve. The bottom panel shows the inferred
radial velocity dispersion with error-bars. The measurements are provided in
Table 6 and Table 7. We observe that the radial velocity very gradually
changes along 200 degrees of the stream, reaching peak values at
$\phi_{1}\approx-60$ deg and $90$ deg. We trace the stream well down to low
Galactic latitudes of $|\phi_{1}|\sim 10$ deg ($\phi_{1}=0$ deg corresponds to
the point where the stream crosses the Galactic plane). Our measurements of
the velocity dispersion in the stream show a typical dispersion of $\sim$ 5
$\mathrm{\,km\,s^{-1}}$ without a clear gradient along the stream. This is
consistent with the measurement by K19 based on SDSS data alone (see their
Figure 10).

###  3.6 Proper motion modelling

We now proceed to extract the stream proper motions by fitting them with a
mixture model similar to the one used for radial velocities. We select stars
by combining previous selections. We use stars from the spectroscopic sample
defined in Section 3.2 that are within 30 km s-1 of the radial velocity track
from Section 3.5, stars within the CMD mask defined in Section 3.4, and stars
that are within 3 degrees of the fiducial stream track on the sky from Figure
1. We then model the distribution of proper motions in a narrow region around
the fiducial proper motion track of the stream. The width of this region is
0.75$\mathrm{\,mas\,yr^{-1}}$ at $|\phi_{1}|<20$ deg and 1.5
$\mathrm{\,mas\,yr^{-1}}$ elsewhere. The reason for modelling a narrow range
of proper motions is that the observed ${S}^{5}$ OC stream candidates were
selected by proper motion (sometimes in a narrow region around the expected
proper motion; see Li et al. 2019 for details); thus to model the full range
of the proper motion distribution we would need to take into account selection
effects, which we avoid here.

The likelihood function for proper motion is similar to previous models and is
a mixture of the stream and background:

$\displaystyle
P(\mu|\phi_{1})=(1-f(\phi_{1}))U(\mu|\mu_{l}(\phi_{1}),\mu_{h}(\phi_{1}))+$
$\displaystyle
f(\phi_{1})N^{*}(\mu|{\mathcal{M}}(\phi_{1}),\mathcal{S_{\mu}}(\phi_{1}))$ (3)

Here, $\mu_{l}$ and $\mu_{t}$ are the $\phi_{1}$-dependent boundaries of the
proper motion fits. $N^{*}$ is the normal distribution truncated at $\mu_{l}$,
$\mu_{t}$, and $f(\phi_{1})$, ${\mathcal{M}}(\phi_{1})$,
${\mathcal{S}_{\mu}}(\phi_{1})$ are the splines for the mixing fraction
between the stream and background, the mean proper motion, and the spread in
proper motions respectively. This model is fit separately to proper motions in
$\phi_{1}$ and $\phi_{2}$ and takes into account the individual proper motion
uncertainties, obtained by projecting the Gaia proper motion covariance
matrices to the $\phi_{1}$, $\phi_{2}$ coordinate system.

The results of the fits are shown in Figure 5. The plot mimics the two bottom
panels of Figure 2 or Figure 4. The top panel of Figure 5 shows the
$\mu_{\phi,1}$ proper motion of individual stars that were modelled with the
grey bands showing the range of proper motions fitted. Red curves show the
samples from the posterior. The second panel from the top shows the inferred
intrinsic spread (standard deviation) of the proper motions in $\phi_{1}$. The
bottom two panels show the same measurements but for $\mu_{\phi,2}$. The
central panel shows the proper motion dispersion along the stream in
$\mathrm{\,km\,s^{-1}}$. The measurements of proper motions and proper motion
dispersions are also given in Tables 8 and 9. The extracted proper motion
tracks are broadly similar to the ones presented in K19, albeit significantly
more accurate. We also notice that a significant intrinsic (i.e. after taking
into account observational errors) proper motion spread in $\mu_{\phi,1}$ is
measured. In the region of $-25\lesssim\phi_{1}\lesssim 25$ deg, where the
stream is at the closest distance to the Sun, the proper motion spread in
$\phi_{1}$ is $\sim 0.15$$\mathrm{\,mas\,yr^{-1}}$: which at 15 kpc
corresponds to 10 $\mathrm{\,km\,s^{-1}}$. That is higher than the dispersion
in the radial velocity of 5 $\mathrm{\,km\,s^{-1}}$ (see Sec 3.5) and is also
higher than the measured velocity dispersion in the $\phi_{2}$ direction of
$\sim 0.05-0.1$$\mathrm{\,mas\,yr^{-1}}$. The difference in measured velocity
dispersions along different directions is somewhat surprising but is in fact
explained by the models; we return to this in Section 6.6.

###  3.7 Stream density model

Figure 6: The spatial density model of the OC stream. The plots represent 2-D
slices of the 4-D model describing the joint distribution of positions and
proper motions. Top panel: The density of stars selected by colour-magnitude
and having proper motions within 0.2 mas yr-1 of the best-fit proper motion
track. Second panel: Best fit stream density on the sky model corresponding to
the same selection of stars as in the top panel. Third panel: The residuals
between the data and the model. Fourth panel: The density of stars in the data
with the best-fit background model subtracted. Bottom panel: The best-fit
model of the stream without the background component of density. In all panels
we only show the sky pixels that were fitted, so pixels with $|b|<5$ and
$E(B-V)>0.25$ are omitted. In all the panels we use linear grey scale, and
with the exception of the middle panel white colour corresponds to zero pixel
values.

As a final step in tracing the OC stream, we proceed to model the stream
density and track across the whole stream extent. Contrasting with the models
described in previous sections, here we work with the entire Gaia EDR3 sample,
rather than focusing on spectroscopic members. We select stars that have
$G<19.5$, small parallaxes ($\omega<3\sigma_{\omega}+0.117$ mas), and that are
within the distance corrected CMD mask from Section 3.4 after shifting them by
the distance modulus from the distance model of Section 3.3. We consider stars
with $-90\,\,{\rm deg}<\phi_{1}<120$ deg, $-7\,\,{\rm deg}<\phi_{2}<9$ deg and
within 0.5 $\mathrm{\,mas\,yr^{-1}}$ from the proper motion track determined
from spectroscopic members.

We note that in this section, we re-model the proper motion track of the
stream. While we could rely on the proper motion track from Section 3.6 and
only model the objects selected to lie close to that track, that would prevent
us from properly modelling the background contamination. Thus we build the
following model $\rho(\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})$ for the
stellar density in bins of $\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2}$. This
is represented by a sum of the background and stream components:

$\displaystyle\rho(\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})=R_{str}(\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})+$
$\displaystyle R_{bg}(\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})$ (4)

The stream density is modelled by a Gaussian cross-section in $\phi_{2}$, and
a Gaussian distribution in proper motions (we note that in contrast to
previous models, we do not use individual proper motion error-bars in the
model: i.e. we model the proper motion distribution as convolved with
observational errors). We describe the mean track on the sky, the logairithm
of the stream width, the logarithm of surface brightness, the proper motion
track, and logarithm of its width by natural cubic splines parametrised by
values at knots:

$\displaystyle
R_{str}(\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})=I(\phi_{1})N(\phi_{2}|\Phi(\phi_{1}),\Sigma(\phi_{1}))$
$\displaystyle
N(\mu_{\phi,1}|{\mathcal{M}}_{1}(\phi_{1})\Sigma_{\mu,1}(\phi_{1}))N(\mu_{\phi,2}|{\mathcal{M}}_{2}(\phi_{1})\Sigma_{\mu,2}(\phi_{1}))$
(5)

Here ${\mathcal{M}}_{1}(\phi_{1})$, ${\mathcal{M}}_{2}(\phi_{1})$,
$\Sigma_{\mu,1}(\phi_{1})$, $\Sigma_{\mu,2}(\phi_{1})$, $I(\phi_{1})$,
${\Phi}(\phi_{1})$, $\Sigma(\phi_{1})$ are the functions for proper motion
tracks ($\mu_{\phi,1}$ and $\mu_{\phi,2}$), width of proper motion tracks,
surface brightness of the stream (per bin), track on the sky, and width of
stream on the sky.

The background is modelled similarly to Erkal et al. (2017) and K19, except
that we add extra terms for the density dependence on proper motions:

$\displaystyle\log R_{bg}($
$\displaystyle\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2})=$ $\displaystyle
B_{000}(\phi_{1})+B_{001}(\phi_{1})\phi_{2}+B_{002}(\phi_{1})\phi_{2}^{2}+$
$\displaystyle\left(B_{100}(\phi_{1})+B_{101}(\phi_{1})\phi_{2}+B_{102}(\phi_{1})\phi_{2}^{2}\right)\mu_{\phi,1}+$
$\displaystyle\left(B_{010}(\phi_{1})+B_{011}(\phi_{1})\phi_{2}+B_{012}(\phi_{1})\phi_{2}^{2}\right)\mu_{\phi,2}$
(6)

We remark that while for simplicity we write Equations 5 and 6 in terms of
$\mu_{\phi,1}$, and $\mu_{\phi,2}$, in practice we model proper motion
residuals with respect to the fiducial proper motion tracks
$M_{1}^{0}(\phi_{1})$, $M_{2}^{0}(\phi_{1})$ from Section 3.2:
$\mu_{\phi,1}^{*}=\mu_{\phi,1}-M_{1}^{0}(\phi_{1})$,
$\mu_{\phi,2}^{*}=\mu_{\phi,2}-M_{2}^{0}(\phi_{1})$.

The model given in Equation 4 is fit to the binned histogram of stars in
$\phi_{1},\phi_{2},\mu_{\phi,1},\mu_{\phi,2}$. We use 2 degree and 0.5 degree
bins for $\phi_{1}$ and $\phi_{2}$ respectively, and $0.033$
$\mathrm{\,mas\,yr^{-1}}$ bins for proper motions. We also excluded pixels
that are very close to the Galactic plane, $-5\,\,{\rm deg}<b<5$ deg and
pixels with high extinction, $E(B-V)>0.25$ from our modelling.

The priors on the model parameters are mostly non-informative and are
available as a part of the Stan model in supplementary materials. The only
weakly informative priors are the priors on the proper motion offsets with
respect to the fiducial proper motion track, ${\mathcal{N}}(0,1)$, the prior
on the logarithm of widths of proper motion track,
${\mathcal{N}}(\log(0.1),0.5)$, and the prior on the logarithm of stream width
on the sky, ${\mathcal{N}}(\log(0.9),0.5)$ (this prior has 16/84 percentiles
of 0.25/2.5 deg, thus covering well the expected range of widths of the OC
stream).

The resulting model has around 300 parameters and is slow to evaluate;
however, as modern versions of Stan allow parallelisation across multiple
cores with the $\tt reduce\\_sum$ functionality, we are still able to sample
the model posterior successfully.

Figure 6 shows the results from the spatial part of the model. Different
panels show the spatial distribution of a subset of fitted stars that lie
close to the proper motion track (top panel), the best fit model (second from
the top panel), the residuals from the model (middle panel), the spatial
distribution of stars with the best-fit background model subtracted (second
from the bottom panel), and the best fit model without the background (bottom
panel). The measurements of the stream track, width, and surface brightness
are provided in Tables 10, 11, and 12. We analyse the details of the linear
density along the stream in Section 4.1.

We see that the model reproduces the data extremely well. The stream is
strongly curved on the sky, both in the North and the South near the stream
edges. We also note significant stream width changes, with the stream being
narrowest at the point of the closest approach to the Sun (and the Galactic
Centre) at $\phi_{1}\sim 10-20$ deg with a width of $\sim$ 0.5 deg. The stream
is noticeably wider in the North at $\phi_{1}\sim 100$ deg, and South at
$\phi_{1}\sim-60$ deg, with widths of $\gtrsim 1$ deg. We note that this
implies quite a drastic difference in the width in physical units (200 pc vs 1
kpc; see Figure 7) since the edges of the stream are significantly further
away at a distance of 60-80 kpc. This on-sky broadening is not expected in a
spherical potential but could arise due to a flattened potential (Erkal et
al., 2016a). Since the width of the stream is related to the progenitor’s
mass, it could also suggest that the progenitor of the OC stream was
significantly more massive when material near the edges of the stream was
stripped than when the material near $\phi_{1}\sim 0$ deg was stripped. We
also note that the measurements of the stream where it crosses the plane are
poorly constrained as the stream there is lost in the contamination. The
spectroscopic sample is somewhat better suited to trace the stream in that
region due to the availability of radial velocities.

Figure 7: OC stream width in kpc as determined from the best fit model shown
in Figure 6. The grey band shows the part of the stream hidden behind the MW
disc, where the width is not directly constrained.

While the model presented in this section gives a good description of stream
density, we point out that K19 modelled the track and density using deeper
DECaLS data, covering part of the Northern sky. That model showed a stream gap
at $\phi_{1}=76$ deg. This gap is not clearly visible in our data and model
because we rely on significantly shallower Gaia data and use 2 degree wide
bins in $\phi_{1}$ which precludes us from seeing small-scale density
features.

Figure 8: The background subtracted histogram of proper motion distribution of
OC stream stars as a function of $\phi_{1}$ based on the Gaia-only data. Here,
we used stars within 1 degree of the stream track on the sky from Figure 6 to
select stream stars, and stars between 2 and 4 degrees away from the stream
track as background stars. The top panel shows $\mu_{\phi,1}$ vs $\phi_{1}$
while the bottom shows $\mu_{\phi,2}$. The curves show the best-fit models
(shifted down by 0.5 $\mathrm{\,mas\,yr^{-1}}$ for clarity). The red curve is
the model fitted to the spectroscopic dataset, and the blue curve is the model
fitted to the Gaia-only data.

As a final check of the proper motion/density model, we verify that the proper
motion track extracted in this section is accurate and matches the one
extracted from spectroscopic data by investigating the proper motion
distribution of likely members from the Gaia-only sample used in this Section.
Figure 8 shows the data-driven proper motion distribution in the stream
extracted from the Gaia sample. To make the figure we take stars within 1
degree of the best-fit stream track as stream stars, and stars that are
between 2 and 4 degrees away from the track as background stars. We then
subtract the appropriately scaled $\phi_{1}$ versus $\mu$ proper motion
distribution of background stars from the distribution of likely stream stars.
That should give a model-independent proper motion stream track. The stream is
clearly visible on both panels as an overdensity. On Figure 8, we also
overplot our best-fit proper motion tracks (shifted for clarity by
0.5$\mathrm{\,mas\,yr^{-1}}$) from the model based on the spectroscopic sample
and the Gaia fit presented in this Section. We see that the models agree well
with each other, giving us confidence in our analysis. The fact that the
models also clearly go over the background subtracted Gaia-only proper-motion
tracks also gives us confidence that the ${S}^{5}$ spectroscopic target
selection based on proper motions did not bias our measurements.

##  4 Analysis of stream measurements

The models fitted to the stream data described in the previous section provide
us with a wealth of information on the full 6-D trajectory of the stream,
stream members, as well as the density of stars along the stream. In this
section, we look at these measurements in more detail.

###  4.1 Stream density

Figure 9: The density of the OC stream as a function of angle along the stream
$\phi_{1}$. The top panel shows the density in solar masses per unit time,
i.e. corrected for the distant-dependent sample incompleteness and for the
accordion effect caused by the faster/slower motion in the
pericentre/apocentre. The bottom panel shows the density only corrected for
the sample completeness. The greyed-out region is the region where the OC
stream crosses the Galactic plane, and stream densities are not well
constrained.

The first property we look at is the stellar density of the stream. While the
map presented in Figure 6 shows that the stream seems ‘brighter’ at some
locations than others, one has to be cautious of several effects that prevent
direct interpretation of this. First, the colour-magnitude selection of stars
used for the modelling is not the same across the stream: the edges of the
stream are significantly further from the Sun and therefore we are only
probing the tip of the red giant branch. To correct for this distance
dependent incompleteness we can assume that the stream consists of a single
stellar population occupying an isochrone with a 12.5 Gyr age,
$\mathrm{\,[Fe/H]}=-2$ and Chabrier IMF (Chabrier, 2005). We use MIST 1.2
isochrones (Choi et al., 2016; Dotter, 2016) and assume that Gaia data is
complete to $G<19.5$ (Boubert & Everall, 2020). Since we have measurements of
the distance to the stream (see Section 3.3) at each $\phi_{1}$ we can use it
correct the stellar density in number counts per degree into stellar mass per
degree.

On top of distance dependent incompleteness, there are also purely kinematic
effects that affect the observed stream density, as parts of the stream that
are closer to us move significantly faster than the stars at the apocentre.
This will lead to an apparent density decrease near the pericentre and a
density increase near the apocentre (i.e. concertina/accordion effect). If we
are interested in the physical stream density changes, for example associated
with the current or previous location of the progenitor along the stream, we
might want to correct for this kinematic stretching/squeezing.

To do this correction, instead of looking at the density of stars per unit
angle $\frac{dN}{d\phi_{1}}$ along the stream, we need to look at the rate
with which stream stars cross a fixed surface perpendicular to the stream,
i.e. $\frac{dN}{d\phi_{1}}\frac{d\phi_{1}}{dt}$555Here the derivative of
$\phi_{1}$ needs to be the reflex corrected proper motion along $\phi_{1}$. .
This essentially describes the rate of flow of stars at each point along the
stream. This is more reflective of the mass loss history of the progenitor
than the stellar density since the stream gets stretched (compressed) as it
speeds up (slows down) along its orbit. In Figure 9, we show two measures of
the density. The top panel shows the linear density of the stream corrected
for the distance-dependent sample incompleteness and accordion effect, while
the bottom panel shows the density corrected just for the distance-dependent
incompleteness.

We notice in the lower panel of Figure 9 that the densities without kinematic
corrections change very significantly along the stream: they are significantly
higher in the South compared to the North. In addition, the density rises by a
factor of almost 10 in the North near $\phi_{1}\sim 100$ deg, which could have
been interpreted as the location of the progenitor. However, all those
features largely disappear in the top panel after correcting for the accordion
effect and where the stream density appears much more uniform. There is some
increase in stream brightness near $\phi_{1}=100$ deg, and possibly at
$\phi_{1}=-60$ deg, but overall the density is surprisingly uniform. This
uniform density is also somewhat disconcerting as it does not help us
constrain where the progenitor of the OC stream is or might have been.

The completeness-corrected linear density map shown in Figure 9 also allows us
to estimate the total stellar mass by integrating along $\phi_{1}$. This gives
us $\sim 5.6_{-0.3}^{+0.6}\times 10^{5}\mathrm{\,M_{\odot}}$ as the mass of
the visible part of the stream (this is well in agreement with an estimate by
Mendelsohn et al. (2022)) . This is equivalent to a luminosity of
$M_{V}\sim-8.2$ assuming $\mathrm{\,[Fe/H]}=-2$ and age = 12.5 Gyr. This
luminosity is somewhat lower (2 $\sigma$) than the luminosity estimate from
K19 based on RR Lyrae count. The luminosity is however in good agreement with
the mean metallicity of $\mathrm{\,[Fe/H]}=-1.9$ (see the next section) based
on the mass metallicity relationship (Kirby et al., 2013), suggesting that we
are not observing a small fraction of the actual stream, or missing a massive
progenitor. The luminosity of the system together with its metallicity make
the progenitor similar in properties to other satellites (classical dwarfs) of
the Milky Way like Sextans and Draco (a similar conclusion was reached based
on a stream model by Hendel et al., 2018).

###  4.2 Metallicity gradients

Given the lack of evidence about the progenitor location based on the stream
density, one may also try to find the location of the stream progenitor
through metallicity gradients. It is a well-known fact that, for example, the
Sagittarius stream shows strong metallicity gradients (Hyde et al., 2015;
Hayes et al., 2020). In fact, it is expected that streams from dwarf galaxies
should show metallicity gradients since most galaxies show negative
metallicity gradients with radius (i.e. more metal-poor outskirts relative to
their centre: Harbeck et al. 2001; Koleva et al. 2011; Okamoto et al. 2017;
Mercado et al. 2021). Since the galaxy outskirts should be stripped first, we
expect the metallicity in progenitor-less streams to show an inverted V shape
as a function of angle along the stream. In other words, the metallicity in
the stream at the location of the progenitor should be highest, and then
decrease moving away from that point (see e.g. Figure 7 in Hayes et al. 2020).
This motivates us to look for possible metallicity gradients in the OC stream.
Here we rely on the catalogue selected based on proper motions, radial
velocities and positions from the main ${S}^{5}$ catalogue. We use the best-
fitting splines from previous sections and also apply additional quality cuts
requiring that the best S/N for each star in both blue and red arms of the
spectra must be >7 (to ensure uniform quality of metallicities).

Figure 10: The metallicities of likely OC stars measured by ${S}^{5}$ vs angle
along the stream. Only stars with high S/N spectra were included. The likely
stream stars were selected based on their proper motion, position on the sky
and radial velocity (black points). The grey points show stars selected
similarly with the exception of a broader proper-motion selection and the
requirement that the velocities differ from the $RV(\phi_{1})$ model by more
than 15 and less than 25 $\mathrm{\,km\,s^{-1}}$. The metallicity and
metallicity gradients shown in the top left/top right of the figure show
measurements from fits to the Southern ($\phi 1<0$)/ Northern($\phi_{1}>0$)
parts of the stream respectively.

We plot metallicity versus $\phi_{1}$ for stars in our spectroscopic sample in
Figure 10. The black points are very likely OC members, while grey points show
the stars selected with the same selection criteria as the members sample, but
with velocity deviating from the stream track of OC stream by more than
15$\mathrm{\,km\,s^{-1}}$ and less than 25$\mathrm{\,km\,s^{-1}}$. This sample
is supposed to represent possible background contamination. Visually, the
metallicities appear roughly constant throughout the stream. To verify that,
we model the whole stream, and the Northern ($\phi_{1}>0$ deg) and Southern
parts ($\phi_{1}<0$ deg) separately, with a Gaussian metallicity distribution
where the mean can linearly change with $\phi_{1}$,

$\mathrm{\,[Fe/H]}=\mathrm{\,[Fe/H]}_{0}+\frac{d\mathrm{\,[Fe/H]}}{d\phi_{1}}(\phi_{1}-\phi_{1,\mathrm{ref}})$

The results of the modelling for the Northern and Southern parts are shown in
the Figure. The modelling of the whole dataset gives an
$\mathrm{\,[Fe/H]}_{0}=-1.9\pm 0.04$, $\mathrm{\,[Fe/H]}$ spread of $0.3$ dex,
and gradient of $\frac{d\mathrm{\,[Fe/H]}}{d100\phi_{1}}=0.0\pm 0.1$ dex/deg.
All the observed gradients are consistent with zero at the $\sim$ 1-sigma
level. This is unfortunate, as similarly to the density along the stream, it
does not shed light on the progenitor location. There are several possible
explanations for this. One is that the progenitor’s location is outside the
range $-80\,\,{\rm deg}<\phi_{1}<100$ deg and therefore we are only seeing the
leading or trailing tail of the stream. Another explanation is that the stream
progenitor may have had an intrinsically small metallicity gradient.

###  4.3 Conservation of quantities along the stream

Figure 11: Changes in the energy and the z-component of the angular momentum
along the OC stream. Each point shows the inferred energy and angular momentum
in 10-degree long segments. The ellipses show the covariances. The colour
indicates the $\phi_{1}$ of the segment; segments with $|\phi_{1}|\leq 20$ deg
were excluded. When computing the covariances, the full 6-D phase space
uncertainties were taken into account, i.e. on the stream track on the sky,
proper motion, distance and radial velocity.

Since the first studies of the MW halo (Helmi et al., 1999), but especially
since the arrival of Gaia data, it became useful to look at halo substructures
in the space of orbital invariants – such as energy, angular momentum and
actions (Myeong et al., 2018; Koppelman et al., 2018; Helmi, 2020) – as
substructures are expected to be more compact there compared to pure phase
space. For example, Bonaca et al. (2021) and Li et al. (2022) looked at a
large group of halo-structures/streams in energy and angular momentum space in
order to associate streams with different progenitors. Also, the lumpiness of
the distribution of stars in the space of conserved quantities was proposed as
a way of constraining the Milky Way potential (Peñarrubia et al., 2012;
Sanderson et al., 2015; Reino et al., 2021).

The energy ($E$) and the z-component of angular momentum ($L_{z}$) are generic
orbital invariants as they are expected to be conserved in time-independent
and axisymmetric potentials. In the case of the OC stream, since we do possess
the full 6-D phase space measurements spanning continuously more than 120 kpc,
we can directly test the conservation of these quantities across a large part
of the Galaxy. Indeed, when Li et al. (2022) divided the OC stream into two
components ($\phi_{1}>0$ deg and $\phi_{1}<0$ deg) in $E-L_{z}$ space, they
found that these two components did not line up. Here, we take a more detailed
look at the change in energy and angular momentum along the stream.

To compute $E$ and $L_{z}$ as a function of $\phi_{1}$, we take stream tracks
of the stream in 6-D: $D(\phi_{1})$, $\phi_{2}(\phi_{1})$,
$\mu_{\phi,1}(\phi_{1})$, $\mu_{\phi,2}(\phi_{1})$, $V(\phi_{1})$ together
with their associated uncertainties measured in Section 3. We then sample from
the uncertainties and evaluate $E$ and $L_{z}$ of stream segments on a grid of
$\phi_{1}$. Here, we use the Milky Way potential from McMillan (2017). Figure
11 shows the resulting energy and $z$-component of angular momentum
measurements at positions along the stream spaced by 10 degrees. The ellipse
shapes show uncertainty/covariance in $E$ and $L_{z}$ (corresponding to
1-$\sigma$ errors). The colour of the ellipse indicates the $\phi_{1}$ angle
of the corresponding stream segment. We note that the uncertainties in $E$ and
$L_{z}$ are often strongly correlated.

Before examining the figure, we also remark that the stream stars will not
necessarily occupy a small region in energy-angular momentum space, as stars
in tidal streams naturally produce a distinctly elongated (bow-tie shaped) E/L
distribution (Yoon et al., 2011; Gibbons et al., 2014). However, in the case
of the OC stream, we expect the intrinsic energy spread in the stream to be of
the order of $v_{\rm peri}\sigma_{v}\approx 0.02\times 10^{5}({\rm km}\ {\rm
s}^{-1})^{2}$ (where $v_{\rm peri}$ is the velocity at the pericentre, $\sim
400$$\mathrm{\,km\,s^{-1}}$, and $\sigma_{v}$ is the velocity dispersion of
the dwarf, $\sim 5$$\mathrm{\,km\,s^{-1}}$). Figure 11 instead shows a very
large variation in both $E$ and $L_{z}$. In particular, $L_{z}$ changes almost
by a factor of two across the stream. The energy variation from the middle of
the stream toward the edges is $0.3\times 10^{5}({\rm km}\ {\rm s}^{-1})^{2}$.
These large changes are clearly not expected from a normal disruption process,
but must be related to the interaction of the stream with the Magellanic
Clouds that causes a spreading of stream stars in energy and angular momentum
space. We look at the model prediction of that distribution in Section 6.6.

###  4.4 Stream direction of motion

Figure 12: The OC stream proper motions corrected for solar reflex motion
compared to the stream track. The grey lines delineate the OC stream as
measured in the paper, while the arrows show the evaluation of the best-fit
proper-motion splines on a grid of $\phi_{1}$. We also show the location of
the LMC to illustrate that the proper motion deviation is towards it. Figure
13: Comparison of stream track orientation and the motion of stars. Top panel:
The stream track direction on the sky $\frac{d\phi_{2}}{d\phi_{1}}$ is shown
in red, while the direction of motion of stars
$\frac{\mu_{\phi,2}\cos\phi_{2}}{\mu_{\phi,1}}$ is shown in blue. The shaded
areas show the 1-sigma uncertainties from sampling the proper motion, track
and distance splines. Bottom panel:  Comparison of the stream track and
stellar velocities in the radial direction. The red line shows the distance
gradient of the stream track along $\phi_{1}$, while the blue line shows the
ratio of radial velocity to the proper motion in $\phi_{1}$ direction. Shaded
bands show the 1-sigma intervals similarly to the top panel. Both panels show
a clear misalignment of the stellar motion and stream direction, meaning that
the stars are not moving along the stream, confirming the effects of the
perturbation by the LMC.

In K19 and E19, it was found that the stars in the OC stream are not actually
moving along the stream track. This is caused by the gravitational pull of the
Large Magellanic Cloud, an effect that has since been detected in several
other streams (Shipp et al., 2021). Here, equipped with the full 6-D track, we
can constrain this behaviour much more effectively. To do so, we use the
proper motion and distance splines inferred in previous sections. Figure 12
shows the solar reflex corrected proper motion vectors along the OC stream
together with the stream track. We notice that, as seen previously, the
velocity vectors are aligned with the stream track for the Northern part of
the stream ($\phi_{1}>0$ deg), while in the South the stream stars are moving
out of the stream plane, pulled towards the LMC. This motion across the stream
is also the reason for the bent-up shape seen in Figure 6, as the Southern
part of the stream is moving up pulled by the LMC.

While Figure 12 provides an illustration of the misalignment of the stars’
motions with that of the stream, we can also directly compare the direction of
motion $\frac{\mu_{\phi,2}\cos\phi_{2}}{\mu_{\phi,1}}$ of the stars versus
stream track direction $\frac{d\phi_{2}}{d\phi_{1}}$. This is shown in the top
panel of Figure 13 (see Fig. 1 of E19, for comparison). Here, we see that the
stream track directions (shown in red) are aligned with the stellar direction
of motion (in blue) for the Northern part of the stream at $\phi_{1}>30$ deg,
with deviations starting to show at lower $\phi_{1}$. On the bottom panel of
Figure 13, we show a similar analysis applied to radial velocity measurements,
i.e. the comparison between the distance gradient in the stream
$\frac{dD}{d\phi_{1}}$ and the radial velocity to proper motion ratio
$\frac{V_{r}\cos\phi_{2}}{\mu_{1}}$. If the stream stars were moving along the
stream track, then those two quantities would be identical; however, the
figure similarly shows that for $\phi_{1}<30$ deg there is a misalignment
between the stars’ motions and the stream track, with the stream distances
decreasing faster than expected from radial velocities alone.

##  5 Dynamical modelling

###  5.1 Stream generation technique and potential models

In order to fit the exquisite data set presented in Section 3 with a dynamical
model, we use a similar technique to E19. In particular, we use the modified
Lagrange Cloud stripping technique from Gibbons et al. (2014) which has been
modified to include the influence of the LMC in E19. Given the significantly
improved data set in this work compared to K19 which E19 fit, we use a Milky
Way and LMC potential with more flexibility. For the baryonic components of
the Milky Way potential, we assume a Hernquist profile (Hernquist, 1990) for
the bulge with a mass of $5\times 10^{9}\mathrm{\,M_{\odot}}$ and a scale
radius of 500 pc, a Miyamoto-Nagai profile (Miyamoto & Nagai, 1975) for the
disc with a mass of $6.8\times 10^{10}\mathrm{\,M_{\odot}}$, a scale radius of
$3$ kpc, and a scale height of $0.28$ kpc based on the MWPotential2014 model
in Bovy (2015). We model the Milky Way’s dark matter halo as an axisymmetric,
generalized NFW (Navarro et al., 1996) with

$\displaystyle\rho(m^{\prime})=\frac{\rho_{0}}{\Big{(}\frac{m^{\prime}}{r_{s}}\Big{)}^{\gamma}\Big{(}1+\frac{m^{\prime}}{r_{s}}\Big{)}^{\beta-\gamma}}\exp\Bigg{(}-\Big{(}\frac{m^{\prime}}{r_{\rm
cut}}\Big{)}^{2}\Bigg{)}\,\,,$ (7)

where

$\displaystyle\rho_{0}=\frac{M_{\rm NFW}}{4\pi r_{s}^{3}}\frac{1}{\ln(1+c_{\rm
NFW})-\frac{c_{\rm NFW}}{1+c_{\rm NFW}}},$ (8)

$M_{\rm NFW}$ is the halo mass, $c_{\rm NFW}$ is fixed to 15 to avoid
degeneracy with the halo mass, $r_{s}$ is the scale radius,
$m^{\prime}=\sqrt{x^{\prime 2}+y^{\prime 2}+z^{\prime 2}/q^{2}}$ is the
flattened radius, $q$ is the flattening of the halo, $\gamma$ and $\beta$ are
the inner and outer slope of the profile respectively, and $r_{\rm cut}$ is
the cutoff radius of the profile. The primed coordinates
($\mathbf{x^{\prime}}$) used here are rotated with respect to the Milky Way
Cartesian coordinates ($\mathbf{x}$),

$\displaystyle\mathbf{x^{\prime}}=R\mathbf{x},$ (9)

so that the axisymmetric halo can be flattened in an arbitrary direction. To
parameterise this rotation, we specify the direction in which the halo is
flattened/stretched in: $(x_{\rm NFW},y_{\rm NFW},z_{\rm NFW})$. This is done
to avoid any period boundaries that would occur if we specified the flattening
direction in terms of polar angles. We convert these coordinates into the
polar angles with

$\displaystyle\theta$ $\displaystyle=\tan^{-1}\left(\frac{y_{\rm NFW}}{x_{\rm
NFW}}\right),$ (10) $\displaystyle\phi$
$\displaystyle=\cos^{-1}\left(\frac{z_{\rm NFW}}{\sqrt{x_{\rm NFW}^{2}+y_{\rm
NFW}^{2}+z_{\rm NFW}^{2}}}\right).$ (11)

With these angles, the matrix $R$ from Equation 9 is specified as

$\displaystyle
R=\begin{pmatrix}\cos\phi\cos\theta&\cos\phi\sin\theta&-\sin\phi\\\
-\sin\theta&\cos\theta&0\\\
\sin\phi\cos\theta&\sin\phi\sin\theta&\cos\phi\end{pmatrix}.$ (12)

We note that an NFW profile would have $\gamma=1$, $\beta=3$, and an infinite
cutoff radius. In this work, we set $r_{\rm cut}=500$ kpc for the MW potential
which is significantly beyond the orbital extent of the OC stream. For
computational efficiency, we evaluate the forcefield of this flattened halo by
rotating into the coordinates where the halo is flattened in the $z$ direction
(i.e. the $\mathbf{x}^{\prime}$ coordinates) and using galpot (Dehnen &
Binney, 1998). Altogether, the Milky Way halo potential has 8 free parameters.

For the LMC, we assume a truncated NFW profile:

$M_{\rm LMC}(<r)=\\\ \begin{cases}M_{\rm LMC}\frac{\log(1+\frac{r}{r_{s,{\rm
L}}})-\frac{r}{r+r_{s,{\rm L}}}}{\log(1+\frac{r_{\rm max,{\rm
L}}}{r_{s}})-\frac{r_{\rm max,\,L}}{r_{\rm max,L}+r_{s,{\rm L}}}},&r\leq
r_{\rm max,\,L}\\\ M_{\rm LMC}\,,&r>r_{\rm max,\,L}\end{cases}$

where $M_{\rm LMC}$ is the total mass of the LMC, $r_{s,{\rm L}}$ is the scale
radius of the LMC, and $r_{\rm max,L}$ is the truncation radius of its NFW
halo. We include this truncation so that we can explore the extent of the
LMC’s dark matter halo. The LMC’s present-day proper motions, radial velocity,
and distance are also free parameters with priors given by observations of
these quantities (Kallivayalil et al., 2013; van der Marel et al., 2002;
Pietrzyński et al., 2019, respectively). Thus, the LMC model has 7 free
parameters.

For the progenitor, we follow E19 and place the progenitor’s present-day
location at $\phi_{1}=6.34$ deg. E19 found that the inferred Milky Way mass
did not depend on the progenitor’s location so we do not vary the progenitor’s
location in our fits. The other parameters describing the present-day phase
space coordinates are left as free parameters: $\mu_{\alpha,\,{\rm
prog}}^{*},\mu_{\delta,\,{\rm prog}},v_{{\rm los},{\rm prog}},d_{\rm
prog},\phi_{2,\,{\rm prog}}$. There are thus 5 free parameters needed to
describe the progenitor. All of our parameters, and their respective priors,
are given in Table 1.

Parameter | Prior | Range
---|---|---
OC progenitor |  |
$\mu_{\alpha,\,\rm prog}^{*}$ | Uniform | ($-6,0$) $\mathrm{\,mas\,yr^{-1}}$
$\mu_{\delta,\,\rm prog}$ | Uniform |  $(0,6)$ $\mathrm{\,mas\,yr^{-1}}$
$v_{r,\,\rm prog}$ | Uniform | ($-250,250$) $\mathrm{\,km\,s^{-1}}$
$d_{\rm prog}$ | Uniform | (0,) kpc
$\phi_{\rm 2,\,prog}$ | Uniform | ($-10$ deg,$10$ deg)
Milky Way |  |
$M_{\rm NFW}$ | Log Uniform | $(1,30)\times 10^{11}\,\mathrm{\,M_{\odot}}$
$r_{s}$ | Uniform |  $(0,)$ kpc
$q_{\rm NFW}$ | Uniform | ($0,1$) or ($1,2$)
$x_{\rm NFW}$ | Normal | $0\pm 1$
$y_{\rm NFW}$ | Normal | $0\pm 1$
$z_{\rm NFW}$ | Normal | $0\pm 1$
$\gamma$ | Uniform | ($0,2$)
$\beta$ | Uniform | ($2,4$)
LMC |  |
$M_{\rm LMC}$ | Log-Uniform | $(10^{8},4\times 10^{11})\,\,\mathrm{\,M_{\odot}}$
$r_{s,\,\rm L}$ | Uniform | (1,) kpc
$r_{\rm max,\,L}$ | Uniform | (8.7,) kpc
$\mu_{\alpha,\,\rm LMC}^{*}$ | Normal |  $1.91\pm 0.02$ $\mathrm{\,mas\,yr^{-1}}$
$\mu_{\delta,\,\rm LMC}$ | Normal |  $0.229\pm 0.047$ $\mathrm{\,mas\,yr^{-1}}$
$v_{r,\,\rm LMC}$ | Normal |  $262.2\pm 3.4$ $\mathrm{\,km\,s^{-1}}$
$d_{\rm LMC}$ | Normal |  $49.59\pm 0.547$ kpc
Noise nuisance parameters |  |
$\sigma_{\phi_{2}}$ | Uniform | ($0^{\circ},1^{\circ}$)
$\sigma_{\rm DM}$ | Uniform | ($0,0.25$)
$\sigma_{\mu_{\phi,1}^{*}}$ | Uniform | ($0,0.25$ $\mathrm{\,mas\,yr^{-1}}$)
$\sigma_{\mu_{\phi,2}}$ | Uniform | ($0,0.25$ $\mathrm{\,mas\,yr^{-1}}$)
$\sigma_{v_{r}}$ | Uniform | ($0,30$ $\mathrm{\,km\,s^{-1}}$)
Table 1: Parameters and corresponding priors for our OC model. We break the
priors down by object. The description of each parameter is provided in
Sections 5.1 and 5.2

###  5.2 Stream generation

Given the progenitor, Milky Way, and LMC parameters, we generate a stellar
stream using the modified Lagrange Cloud stripping (mLCs) method of Gibbons et
al. (2014) which E19 generalized to include an LMC. We include the dynamical
friction from the Milky Way on the LMC using the results of Jethwa et al.
(2016). We also account for the motion of the Milky Way in response to the LMC
(e.g. Gómez et al., 2015) which is known to affect many streams (e.g. Erkal et
al., 2019; Vasiliev et al., 2021; Ji et al., 2021). This is done by modelling
the Milky Way and LMC as individual particles sourcing their respective
potentials.

We model the stream progenitor as a Plummer sphere (Plummer, 1911) with an
initial mass of $2.67\times 10^{7}\mathrm{\,M_{\odot}}$ and a fixed scale
radius of 1 kpc. Since we do not see the progenitor in the data, the mass of
the progenitor is modelled to change linearly from the initial mass to zero at
the present day. These parameters were selected to approximately reproduce the
width of the stream on the sky. Before generating the stream in the
simulation, the progenitor is first rewound for 4 Gyr in the presence of the
Milky Way and LMC, and then the system is evolved forwards. Tracer particles
are released from the progenitor at the Lagrange points using the mLCs method.

###  5.3 Likelihood and MCMC exploration

In order to compare each stream model with the data (i.e. sky track, distance
modulus, proper motions, and radial velocity), we compute the likelihood of
the data given the model. As is the case for our data, we use values at spline
knots and assume the measurements are uncorrelated. For each observable, we
define the likelihood at each knot location as

$\log\mathcal{L}_{i}=-\frac{1}{2}\log\Big{(}2\pi\left(\sigma_{i,\rm
data}^{2}+\sigma_{i,\rm sim}^{2}\right)\Big{)}-\frac{1}{2}\frac{\left(m_{i,\rm
data}-m_{i,\rm sim}\right)^{2}}{\sigma_{i,\rm data}^{2}+\sigma_{i,\rm
sim}^{2}},$

where $m_{i,\rm data}$ is the measurement at the $i^{\rm th}$ knot location,
$\sigma_{i,\rm data}$ is the uncertainty on this value, and $m_{i,\rm sim}$ is
the value in the simulation. To infer the simulated value, we fit a straight
line to the simulated stream particles within 5 deg of the knot location to
infer the mean ($m_{i,\rm sim}$) and uncertainty on the mean ($\sigma_{i,\rm
sim}$). We sum the log-likelihood terms over all knot locations and
observables to obtain the full likelihood function of the data given the
simulation.

We note that our model has shot noise since it is represented with a finite
number of particles. To ensure that the likelihood function is sufficiently
smooth, we require that the shot noise of our stream model is small compared
to the observed uncertainties ($\sigma_{i,\rm sim}<\sigma_{i,\rm data}/5$). In
order to achieve this, we strip 120,000 particles per pericenter. This number
was chosen based on initial testing to meet our uncertainty requirement. We
have also checked with the posterior chains and found that it is met $\sim
99\%$ of the time.

Based on initial testing, we found that our best-fit models were not able to
fully explain the data, and the typical $\chi^{2}$ per degree of freedom of
our best-fit models is $\sim 2.5$. This can be caused by either inadequacy of
our model, or possibly underestimated uncertainties (or ignored correlations
in measurements).

To account for this potential problem, we introduced an extra noise nuisance
parameter ($\sigma$) for each observable, which is added in quadrature with
the errors to give a likelihood of

$\displaystyle\log\mathcal{L}_{i}=-\frac{1}{2}\log\Big{(}2\pi\left(\sigma_{i,\rm
data}^{2}+\sigma_{i,\rm sim}^{2}+\sigma^{2}\right)\Big{)}$
$\displaystyle-\frac{1}{2}\frac{\left(m_{i,\rm data}-m_{i,\rm
sim}\right)^{2}}{\sigma_{i,\rm data}^{2}+\sigma_{i,\rm sim}^{2}+\sigma^{2}}.$
(13)

We add one nuisance parameter for each type of observable, giving us five
additional parameters. The maximum of the prior on each parameter is chosen to
be slightly bigger than the largest residual between the model and data based
on initial testing. The final number of parameters is 25 in the case of the
oblate and prolate halo, and 21 parameters for a spherical halo.

To sample the posterior obtained by combining the priors shown in Table 1 and
the likelihood function described above, we use the Markov Chain Monte Carlo
(MCMC) sampler emcee (Foreman-Mackey et al., 2013). We use 250 walkers for
7,500 steps, with a burn-in of 3,750 steps. The post-burn-in posterior chains
are available through Zenodo (see Data Availability section).

Figure 14: Best-fit OC model compared to the data. The stream measurements are
shown with red/grey error-bars and the stream model is shown with blue points.
The grey points are not included in the modelling since these occur in a
region where the model predicts a kink on the sky, i.e. a large range of
$\phi_{2}$. From top to bottom, the panels show the stream on the sky, in
distance modulus, proper motions, and radial velocity. The model is a good
representation of the data, except for the radial velocity in the right-most
bin. A movie showing the evolution of these observables is available here.

The key parameter measurements from the posterior samples for the three MW
halo models we consider (spherical, oblate, prolate) are given in Table 2. The
oblate halo model has the highest log-likelihood and thus represents our best
model for the stream, as well as the Milky Way and LMC. In the following
discussion, we will take the oblate halo as our fiducial model. For
completeness, we also give the nuisance parameters for these three fits in
Table 13. The best-fit stream model is compared with the data in Figure 14.

Parameter | Spherical | Prolate | Oblate
---|---|---|---
Milky Way |  |  |
$M_{\rm NFW}$ ($10^{11}\mathrm{\,M_{\odot}}$) | $8.41^{+5.11}_{-3.74}$ | $13.78^{+2.88}_{-4.26}$ | $11.22^{+4.77}_{-3.73}$
$r_{s}$ (kpc) | $13.19^{+6.66}_{-4.96}$ | $10.63^{+3.65}_{-2.16}$ | $16.10^{+5.97}_{-5.17}$
$q_{\rm NFW}$ | 1 | $1.40^{+0.12}_{-0.10}$ | $0.55^{+0.10}_{-0.08}$
$l_{\rm NFW}$ | - | $101.49^{+9.48}_{-12.67}$ | $16.19^{+4.83}_{-8.27}$
$b_{\rm NFW}$ | - | $35.99^{+11.26}_{-9.39}$ | $32.98^{+5.19}_{-6.81}$
$\gamma$ | $1.11^{+0.37}_{-0.45}$ | $0.66^{+0.35}_{-0.38}$ | $1.27^{+0.31}_{-0.39}$
$\beta$ | $3.19^{+0.36}_{-0.46}$ | $3.82^{+0.13}_{-0.23}$ | $2.97^{+0.32}_{-0.41}$
$M_{\rm MW}$ (32.4 kpc) | $2.97^{+0.15}_{-0.14}$ | $3.06^{+0.10}_{-0.11}$ | $2.85^{+0.09}_{-0.09}$
$M_{\rm MW}$ (50 kpc) | $3.90^{+0.22}_{-0.23}$ | $3.81^{+0.13}_{-0.13}$ | $3.77^{+0.26}_{-0.19}$
$M_{\rm vir,\,MW}$ ($10^{11}\mathrm{\,M_{\odot}}$) | $7.19^{+3.56}_{-1.12}$ | $5.68^{+0.40}_{-0.26}$ | $7.74^{+4.14}_{-1.48}$
$r_{\rm vir,\,MW}$ (kpc) | $236^{+34}_{-13}$ | $219^{+5}_{-3}$ | $242^{+37}_{-17}$
LMC |  |  |
$M_{\rm LMC}$ ($10^{10}\mathrm{\,M_{\odot}}$) | $10.28^{+3.51}_{-3.47}$ | $11.70^{+2.79}_{-2.60}$ | $12.85^{+2.76}_{-2.34}$
$r_{s,\,\rm L}$ (kpc) | $2.71^{+2.59}_{-1.18}$ | $3.84^{+1.86}_{-1.56}$ | $2.48^{+1.76}_{-1.06}$
$r_{\rm max,\,L}$ (kpc) | $134.83^{+44.31}_{-82.46}$ | $152.40^{+33.10}_{-44.31}$ | $156.22^{+32.56}_{-41.26}$
$M_{\rm LMC}$(32.8 kpc) ($10^{10}\mathrm{\,M_{\odot}}$) | $5.75^{+0.89}_{-0.71}$ | $5.95^{+0.63}_{-0.64}$ | $7.02^{+0.99}_{-0.86}$
$\Delta\log\mathcal{L}$ | $-23.2$ | $-9.8$ | 0
Table 2: Key parameter measurements from the posterior samples. In addition to
the model parameters, we also provide some derived quantities such as masses
within fixed apertures and the virial mass of the MW. We show parameters for
three models of the MW halo: oblate, prolate, and spherical. Note that the
values correspond to the median with uncertainties from the 84.1 and 15.9
percentiles.

###  5.4 Milky Way potential measurement

In this section, we look at the inference of the Milky Way potential from the
posterior samples presented in Section 5.3.

In order to highlight our mass constraints, we first compute the enclosed mass
profile of the Milky Way. This is done by taking individual posterior samples,
computing the resulting enclosed mass profile, and then computing the 15.9,
50, and 84.1 percentiles of the enclosed mass profile at each radius. We note
that due to the flattening of the mass distribution, we numerically integrate
the density profile in order to compute the mass profile. The resulting
enclosed mass profile of the Milky Way is shown in the top panel of Figure 15.
Note that this enclosed mass includes the bulge, disc, and dark matter halo.
For comparison, we show several other recent mass measurements. Our enclosed
mass profile is a good match to the other measurements, especially those which
fit stream data with dynamical stream models such as Küpper et al. (2015);
Erkal et al. (2019); Vasiliev et al. (2021). We also highlight the
Galactocentric distance range spanned by our stream data.

Next, we consider the fractional uncertainty of the enclosed dark matter mass
in order to determine the region in which we best measure the stream data.
This is shown in the bottom panel of Figure 15. Such a fractional uncertainty
was originally proposed by Bonaca & Hogg (2018). Using a suite of mock
streams, they showed that given a suitably flexible potential, the fractional
uncertainty in the acceleration (or equivalently the mass enclosed) should be
minimized at the present-day location of the stream. For the first time, we
are able to show the same result with real data. In particular, we best
measure the Milky Way’s dark matter mass at a radius of 32.4 kpc, with a mass
of $M_{\mathrm{DM}}=(2.18\pm 0.1)\times 10^{11}\mathrm{\,M_{\odot}}$,
corresponding to an uncertainty of 4.2%, and an enclosed mass of
$M_{\mathrm{MW}}(r<32.4{\rm kpc})=(2.85\pm 0.1)\times
10^{11}\,\mathrm{\,M_{\odot}}$. This sits roughly in the middle of the
present-day Galactocentric range spanned by the observed parts of the stream.

Figure 15: Constraints on the enclosed mass profile of the Milky Way. Top
panel shows the enclosed mass (including all the Galactic components) as a
function of radius. The red curve shows the median enclosed value, while the
dark and light red-shaded regions show the 1 and 2$\sigma$ confidence
intervals. The vertical dashed lines show the radial extent of our data set.
The observations are Deason21 (Deason et al., 2021), Erkal19 (E19), Vasiliev21
(Vasiliev et al., 2021), McMillan17 (McMillan, 2017), Gibbons14 (Gibbons et
al., 2014), Kupper15 (Küpper et al., 2015), CorreaMagnus22 (Correa Magnus &
Vasiliev, 2022) and Reino21 (Reino et al., 2021). Bottom panel shows the
fractional uncertainty on the enclosed dark matter mass, which demonstrates
that the mass is best measured in the region where the OC stream is observed
and that the uncertainties are substantially larger outside of this region.
The most precise measurement of the Milky Way’s enclosed mass is at 32.4 kpc,
with $M_{\rm MW}({\rm r}<32.4{\rm kpc})=(2.85\pm 0.1)\times
10^{11}\mathrm{\,M_{\odot}}$, which corresponds to a fractional uncertainty of
4.2% on the enclosed dark matter mass. The constraints on the mass profile
shown on the panel of the figure are provided in digital form in supplementary
materials.

Next, we compare the enclosed dark matter mass profile of our fiducial
(oblate) fits with the prolate and spherical halo fits in Figure 16. In
particular, we show the fractional difference between each enclosed mass
profile and the oblate mass profile. We see that over the radial range spanned
by the data, our different fits all produce a consistent enclosed mass
profile. This shows the robustness of our inferred mass profile. However, we
note that beyond this radial range, the different fits give different results.

As a demonstration of this, we compute the virial mass and radius using the
standard definition from Bryan & Norman (1998): i.e. the mass enclosed by a
radius within which the average density is $\Delta_{c}\rho_{\rm crit}$, where
$\rho_{\rm crit}=3H_{0}^{2}/8\pi G$ is the critical density of the Universe,
$H_{0}=67.37$ km s-1 Mpc-1 is the present-day Hubble constant, and the pre-
factor $\Delta_{c}\sim 103$, with cosmological parameters from Planck
Collaboration et al. (2020). We present the inferred virial mass of the Milky
Way in Table 2. While the oblate and spherical haloes have similar virial
masses, the mass inferred from the prolate halo is significantly smaller. We
stress that the virial masses are extrapolations and that the most robust
measurement is within the radial range where the OC stream is observed.

Figure 16: Comparison of the Milky Way’s enclosed mass for the models of the
MW DM halo. The y-axis shows the fractional difference between the dark matter
mass inferred from each fit and that inferred from the oblate fit. The
enclosed mass of the Milky Way inferred from these fits is consistent at the
$1\sigma$ level over most of the radial range considered, showing that the
inferred mass is robust.

The flexible models used in this work also allow us to measure the density
slope of the inner ($\gamma$) and outer ($\beta$) Milky Way dark matter halo,
as well as the flattening of the Milky Way halo ($q_{\rm NFW}$).
Interestingly, we see that the oblate and spherical halo have inner and outer
slopes consistent with an NFW profile, i.e. consistent with 1 and 3
respectively. In contrast, the prolate halo has a much steeper outer profile,
$\beta_{\rm prolate}=3.82^{+0.11}_{-0.12}$. This rapid fall-off in density is
why the prolate halo has a much smaller virial mass. In addition, we find that
both the oblate $\left(q_{\rm NFW}=0.56^{+0.10}_{-0.08}\right)$ and prolate
$\left(q_{\rm NFW}=1.40^{+0.12}_{-0.10}\right)$ halo are preferred over the
spherical halo. We discuss the oblate and prolate halo flattening in more
detail in Section 6.2.

###  5.5 LMC mass measurement

We now explore the inferred LMC potential. As with the Milky Way potential in
Section 5.4, we consider the enclosed mass profile of the LMC. We note that
previous works which used streams to measure the mass of the LMC (e.g. Erkal
et al., 2019; Shipp et al., 2021; Vasiliev et al., 2021) all considered a
family of single parameter models for the LMC. In particular, Erkal et al.
(2019); Shipp et al. (2021) assumed a Hernquist profile described by the total
mass with a scale radius chosen to satisfy the mass constraint of van der
Marel & Kallivayalil (2014) at 8.7 kpc. Similarly, Vasiliev et al. (2021)
assumed a smoothly truncated NFW profile described by the total mass with a
scale radius chosen to match the mass enclosed in the inner $\sim 9$ kpc of
the LMC (see their Fig. 3).

In comparison, our flexible model does not put priors on the enclosed mass
profile of the LMC. Thus, for the first time, we present a measurement of the
enclosed mass profile of the LMC in the top panel of Figure 17. For reference,
we show the closest approach distance of the OC stream to the LMC from our
best-fit model. Within this radius, we should not be sensitive to the mass
distribution. We also show mass measurements of the LMC based on rotation
curve measurements for comparison.

Motivated by the results of Bonaca & Hogg (2018), we show the fractional
uncertainty of the enclosed LMC mass in the bottom panel of Figure 17. We see
that the LMC mass is best measured where we broadly expect, i.e. beyond the
closest approach (within which it is only sensitive to the total mass), but
still near the closest approach (where the velocity kicks are the largest). In
particular, there is a large uncertainty on the mass within the closest
approach distance, within which we are not sensitive. Similarly, at large
distances, the enclosed mass is poorly measured because in this region the
acceleration from the LMC is small – and smaller than the acceleration from
the Milky Way. Between these two limits, we achieve the best measurement of
the LMC’s enclosed mass at 32.8 kpc of $M_{\rm LMC}(<32.8\,{\rm
kpc})=7.02^{+0.99}_{-0.86}\times 10^{10}\mathrm{\,M_{\odot}}$.

Figure 17: Constraints on the enclosed mass profile of the LMC. Top panel
shows the enclosed mass profile of the LMC compared to observations. The
vertical dashed grey line shows the median closest approach distance of the OC
stream to the LMC (21.3 kpc) for our models, restricted to the same $\phi_{1}$
range as the data, within which we do not expect to be sensitive. The previous
measurements for the inner LMC are van der Marel14 (van der Marel &
Kallivayalil, 2014), Cullinane20 (Cullinane et al., 2020), and Wan20 (Wan et
al., 2020). The previous measurements for the total LMC mass are Erkal19
(E19), Vasiliev21 (Vasiliev et al., 2021), and Shipp21 (Shipp et al., 2021).
Note that these outer measurements are staggered radially for readability.
Bottom panel shows the fractional uncertainty in the mass. The most precise
measurement of the LMC’s enclosed mass is made at 32.8 kpc with $M_{\rm
LMC}(<32.8\,{\rm kpc})=7.02\times 10^{10}\mathrm{\,M_{\odot}}$, which
corresponds to a fractional uncertainty of 13.2%. The constraints on the LMC
mass profile shown on the top panel of the figure are provided in digital form
in supplementary materials.

Interestingly, Figure 17 also shows that our fits disfavour compact LMC mass
profiles, i.e. where all of the mass sits within the closest approach
distance. Instead, our models prefer an extended dark matter halo for the LMC
where different parts of the stream experience different LMC masses. In order
to better showcase this, we show the inferred density profile of the LMC in
Figure 18. This shows that our models require an LMC whose dark matter halo
extends out to at least $\sim 53$ kpc at the 2$\sigma$ level. Such an extended
dark matter halo is in agreement with the recently observed extent of the LMC
hot corona (Krishnarao et al., 2022) and is consistent with the LMC being on
its first approach to the Milky Way (e.g. Besla et al., 2007).

Figure 18: Inferred density of the LMC. The red curve shows the median density
and the dark (light) shaded regions show the $1\sigma$ (2$\sigma$) confidence
intervals. Interestingly, our fits prefer an extended LMC halo that extends
beyond $\sim 53$ kpc.

##  6 Discussion

###  6.1 Stream dispersions

In Sections 3 & 4, we showed that with the combination of Gaia and the
${S}^{5}$ survey, we can measure the 6-D track and dispersions of the stream.
We also showed in Section 4 that we reliably detect the dispersion of proper
motions along the stream which appears to be around 10
$\mathrm{\,km\,s^{-1}}$: higher than the dispersion of line-of-sight
velocities of 5 $\mathrm{\,km\,s^{-1}}$. In order to highlight the reliability
of the dispersion measurements, as well as to further validate our stream
models, we compare the observed and simulated dispersions in Figure 19. As a
caveat, we note that the mass of the progenitor was selected to reproduce the
average stream width on the sky. That said, the model is able to reproduce the
trends seen in the data, such as the rise and fall of the dispersion in the
proper motion along the stream, as well as all of the other dispersions.
Interestingly, at large $\phi_{1}$, the width of the distance modulus and the
stream on the sky are both larger in the data than in the model. This may be
due to a different mass loss history for the real stream (e.g. it was more
massive in the past than in our model) or some additional complexity in the
potential which could increase these dispersions (e.g. Erkal et al., 2016a).

Figure 19: Width of OC stream in each observable. In each panel, the blue line
shows the dispersion of the best-fit model from Figure 14. From top to bottom,
the panels show the dispersion of the stream on the sky, in distance modulus,
in proper motions, and in radial velocity. In each panel, the red error bars
show the observed dispersions and the blue lines show the dispersion of the
best-fit stream model for comparison. While the model was tuned to match the
mean width on the sky, it is able to match the other dispersions.

###  6.2 Oblate versus prolate halo

In Section 5.4, we discussed the Milky Way potential measurement with a focus
on the mass enclosed, which is robust across our choice of halo flattening. In
this section, we explore why our fits prefer both a prolate and oblate halo
over a spherical halo. First, in Figure 20, we focus on what features in the
data the oblate and prolate haloes are able to represent. We highlight the on-
sky track in this plot since this is the main improvement in the fit; i.e. the
improvement in the likelihood between the spherical and flattened halo fits is
mainly driven by an improved match to the on-sky track (as in E19, ). In
particular, the oblate and prolate halo are able to replicate the downturn in
the on-sky track at $\phi_{1}>100$ deg. In contrast, the spherical halo
provides a much shallower downturn.

Figure 20: Best-fit stream models for different halo models. From top to
bottom, the panels show the best-fit stream in an oblate, prolate, and
spherical halo respectively. Compared to the spherical halo, the main
improvement of using an oblate or prolate halo is that it allows us to
reproduce the observed down-turn in the stream track at $\phi_{1}>100$ deg.

Based on the downturn on the sky, we next investigate how the acceleration
fields of the oblate and prolate halo differ. In order to study this, we
consider the acceleration field from the Milky Way halo in the plane of the
stream. We define the stream plane using the rotation matrix from K19, which
defines the stream coordinates and hence the stream plane which goes through
the Sun. One advantage of this coordinate system is that a force which points
down in these coordinates (i.e. in the negative $z$ direction), will also
point down on the sky in the OC coordinate system (i.e. in the negative
$\phi_{2}$ direction). We shift this plane to go through the centre of the
Galaxy. Given this stream plane, we evaluate the ratio of the out-of-plane
acceleration to the radial acceleration for points on an evenly spaced
Cartesian grid in this plane. We average this acceleration ratio across 100
random draws from the posterior chains and show the results in Figure 21.
Interestingly, the oblate and prolate halo fits have a similar amplitude and
orientation of the acceleration ratio in this plane. Furthermore, in the
region where the stream is bending down on the sky (i.e. $\phi_{1}>100$ deg),
both haloes have an acceleration field that points down in the stream plane.
Thus, despite their seemingly disparate shapes, both the oblate and prolate
halo fits seem to be trying to reproduce the same acceleration field which can
bend the stream down on the sky.

Figure 21: Ratio of out-of-plane to radial acceleration for the oblate and
prolate halo fits. The left and right panels show the ratio of the forces out
of the stream plane to the radial force in the stream plane. This shows that
the oblate and prolate halo fits create a broadly similar force field in the
plane of the OC stream. The black points show the best-fit stream model.

Finally, we compare the orientation of the oblate and prolate haloes in Figure
22. In previous fits to the OC stream, E19 noted that the oblate halo
orientation was consistent with the orientation of the LMC’s orbital angular
momentum vector, and similarly that the prolate halo orientation was
consistent with the present-day vector to the LMC, suggesting a possible
connection between the halo shapes and the LMC. We show the results of E19, as
well as the LMC’s angular momentum vector and present-day position vector for
comparison in Figure 22. We see that while the prolate halo is still aligned
with the LMC’s present-day position, the oblate halo is now slightly
misaligned with the LMC’s angular momentum vector. For comparison, we also
show the orientation of the short-axis of the halo inferred with the
Sagittarius stream by Vasiliev et al. (2021). In that work, they considered a
triaxial halo but the intermediate and long axes were nearly identical,
consistent with an axisymmetric, oblate halo. With the improved accuracy on
the potential in this work, we see that the Milky Way potential shapes
inferred with the OC stream and the Sgr stream are now inconsistent. In light
of this discrepancy, and also the broadly consistent acceleration field of the
oblate and prolate haloes in the OC stream plane, we note that these fits may
be trying to describe the deformations of the Milky Way and LMC haloes which
are known to affect the OC stream (Lilleengen et al., 2023).

Figure 22: Comparison of the orientation of the oblate and prolate halo
inferred in this work. The top (bottom) panel shows the location of the short
(long) axis of the oblate (prolate) halo respectively, as viewed from the
Milky Way centre. The red error bars show the orientation of the oblate and
prolate halo inferred in E19. The blue points show the angular momentum vector
of the LMC (top panel) and the position of the LMC (bottom panel). Note that
the errors in these quantities are smaller than the size of the blue point. In
the top panel, we also show the orientation of the short-axis from Vasiliev et
al. (2021) from a fit to the Sagittarius stream.

###  6.3 Perturbations from classical satellites

In Section 5, we focused on the combined influence of the Milky Way and the
LMC on the OC stream. Next, we assess whether any of the other classical
satellites also could have perturbed the OC stream. In order to do this, we
take the same stream generation machinery as in Section 5.2 and add an
additional perturber in the form of a Plummer sphere with a mass of
$10^{9}\mathrm{\,M_{\odot}}$ and a scale radius of 1 kpc. For the present-day
location of the dwarfs, we use the distances and radial velocities from
McConnachie (2012). For proper motions, we use the results of Gaia
Collaboration et al. (2018) for all the classical dwarf galaxies except Leo I,
Leo II, and the SMC. For these three, we use proper motion measurements from
Sohn et al. (2013), Piatek et al. (2016), and Kallivayalil et al. (2013)
respectively.

For the parameters of the Milky Way, LMC, and OC stream, we take the best-fit
model from the oblate halo fit described in Section 5.3. For each of the
classical dwarfs, we sample its present-day phase-space coordinates 10 times
and rewind it self-consistently in the combined presence of the Milky Way,
LMC, and OC stream progenitor. The system is then evolved forwards while stars
are ejected from the OC progenitor to form the OC stream. We then compare the
resulting OC stream with the one without a perturbation from the classical
satellite (i.e. the fiducial stream). In order to evaluate the changes to the
stream observables, we make mock observations of these streams at the
locations of the knots and compute the difference in the observables with the
fiducial stream. We then compare these differences with the observed
uncertainty at each knot to determine how significant they are. We find that
the only satellite which results in a $>3\sigma$ change in the observables is
the SMC, which perturbs the observed distance by $\sim 3\sigma$ in 2 out of
the 10 realizations. As a result, we find that the OC stream is likely not
strongly affected by the other classical satellites given the present-day
data.

###  6.4 Association with dwarfs and globular clusters

Next, we search for associations between known MW satellites (ultra-faint
dwarfs and GCs) with the OC stream. Previous works have searched for such
associations by comparing the present-day phase space coordinates of each
satellite with observables of the OC stream (e.g. K19, ). In contrast, we
integrate the orbits of each dwarf and GC backwards in the presence of the
stream for 4 Gyr and compute how close each satellite passes to stream
particles in phase space. In particular, for each satellite, we simultaneously
sample from its present-day phase space coordinates and sample the posterior
chains for the OC stream and potential parameters. This sampling is performed
100 times for each satellite.

For the ultra faint dwarfs, we use the sample of 46 ultra faints in Pace et
al. (2022) that have full 6-D phase-space coordinates. In Figure 23, we show
the closest approach distance and relative velocity at this closest approach
for the dwarfs that have an approach within 10 kpc and 100
$\mathrm{\,km\,s^{-1}}$. Interestingly, we see that Grus II stands out as
having realizations that have close passages ($\lesssim 0.1$ kpc) with the OC
stream particles at low relative velocities ($\lesssim 50$
$\mathrm{\,km\,s^{-1}}$). This raises the possibility that Grus II could have
been bound to the progenitor of the OC stream. For reference, we show the
escape velocity curves of a $10^{8}\,\mathrm{\,M_{\odot}}$ and
$10^{9}\,\mathrm{\,M_{\odot}}$ subhalo which are modelled as Hernquist
profiles, with a halo mass-scale radius relation taken from Erkal et al.
(2016b). Interestingly, we find that 11% of the realizations would have been
bound to a $10^{9}\,\mathrm{\,M_{\odot}}$ progenitor.

We note that during each rewinding, Grus II has multiple close approaches to
the OC stream. The most recent of these occur $\sim 300-400$ Myr ago, with
more ancient passages $\sim 1.75-2.25$ Gyr ago. Interestingly, in some of
these earlier passages, Grus II passes very close to the progenitor we have
used in this work. Thus, these rewindings are consistent with the picture
where Grus II could have been bound to the original host of the OC stream
$\sim 2$ Gyr ago, which subsequently disrupted. Thus, Grus II would be the
first ultra-faint dwarf galaxy associated with a classical dwarf galaxy.

We repeat the same analysis for the sample of 160 Milky Way globular clusters
using the 6-D catalogue tabulated by Vasiliev & Baumgardt (2021) and Baumgardt
& Vasiliev (2021). As with the dwarfs, we sample each globular cluster 100
times while simultaneously sampling the stellar stream model. Unlike the
dwarfs, we do not find any close and low-velocity passages: i.e. no globular
cluster passes within 10 kpc and 100 $\mathrm{\,km\,s^{-1}}$ during the
rewinding.

Figure 23: Closest approach distance and relative velocity of each dwarf with
the OC stream. The black solid (dashed) curves show the escape velocity of a
$10^{9}\mathrm{\,M_{\odot}}$ ($10^{8}\mathrm{\,M_{\odot}}$) dark matter halo.
Grus II stands out as having a high probability of having a very close
encounter with the OC stream. In particular, there are many realisations where
Grus II could have been bound to a $10^{9}\mathrm{\,M_{\odot}}$ host halo of
the OC progenitor.

###  6.5 Interaction with LMC

Of the streams modelled to date, the OC stream has one of the closest and
strongest encounters with the LMC (Shipp et al., 2021). In order to better
understand this encounter, as well as what we can learn from it, we show the
closest approach distance and time of the closest approach in Figure 24 for
the best-fit model from Figure 14. In particular, we see that the kink in the
stream track at $\phi_{1}\sim-80$ deg corresponds to the point of closest
approach with the LMC. We note that previous fits to the OC stream also
predicted such a kink (Erkal et al., 2019; Shipp et al., 2021). For our best-
fit model, this approach is very close with some of the particles passing
within $\sim 6$ kpc of the LMC. Furthermore, this passage occurs $\sim 370$
Myr ago, when the LMC was at a distance of $\sim 107$ kpc from the Milky Way.
This suggests that the OC stream should be a sensitive probe of the LMC’s past
location.

In Figure 25, we show how well the fits in this work constrain the past orbit
of the LMC. In order to show the uncertainty in the LMC’s past orbit without
using our OC stream fit, we rewind the LMC in the Milky Way potential from
McMillan (2017). In particular, we sample 100 Milky Way potentials from the
posterior chains in that work. For each of these realizations of the Milky Way
potential, we also sample the LMC’s proper motion, radial velocity, and
distance from their observed values (Kallivayalil et al., 2013; van der Marel
et al., 2002; Pietrzyński et al., 2019, respectively). We model the LMC as a
Hernquist profile (Hernquist, 1990) and sample its mass from the value
measured in E19, $M_{\rm LMC}=(1.38\pm 0.26)\times
10^{11}\mathrm{\,M_{\odot}}$, with a scale radius such that the circular
velocity at 8.7 kpc matches the observed value of 91.7 $\mathrm{\,km\,s^{-1}}$
(van der Marel & Kallivayalil, 2014). We show the orbit for the past 500 Myr
and show the location at a lookback time of 254 Myr, motivated by the median
closest approach of the LMC and the OC stream in our models (see Fig. 24). We
also show the uncertainty of the LMC’s past orbit from the fits in this work
by taking 100 samples from our posterior chains. Figure 25 shows that fits to
the OC stream in this work give a smaller uncertainty in the past trajectory
of the LMC. In addition, we see that the LMC’s orbit in the potential from
McMillan (2017) remains closer to the Milky Way. This is due to the larger
mass of the Milky Way in the models of McMillan (2017, see Fig 15).

In order to further explore how sensitive the OC stream is to the LMC’s past
trajectory, we perform two additional tests. First, we take the best-fit OC
stream model (i.e. from Fig. 14) and turn off the dynamical friction the LMC
experiences in the presence of the Milky Way. The resulting stream is a poor
fit to the data with a difference in the likelihood of $\Delta\log
L\sim-3500$. Without dynamical friction, the LMC’s location during the closest
approach (254 Myr ago) changes by $\sim 6$ kpc, which is comparable to the
LMC’s closest approach distance to the observed regions of the OC stream
($\sim 21.3$ kpc). This shows that current observations are already sensitive
to the dynamical friction that the LMC has experienced.

As a second test, we re-fit the OC stream with an extra parameter that can
scale up or down the dynamical friction the LMC experiences. In particular, we
use the same setup as in Section 5.3 for the oblate Milky Way halo, with one
extra parameter ($\lambda_{\rm DF}$). This parameter multiplies the dynamical
friction term the LMC experiences and it has a log uniform prior between
$10^{-3}$ to 10. As expected from the first test, this parameter is strongly
constrained by the existing data with a posterior of $\log_{10}\lambda_{\rm
DF}=0.31^{+0.20}_{-0.48}$ consistent with $\lambda_{DF}=1$.

Finally, we note that further work is needed to determine whether this is a
robust measurement of dynamical friction. In particular, in the model used for
this paper, we have neglected the deformations of both the Milky Way and LMC
dark matter haloes which have been predicted to affect the OC stream
(Lilleengen et al., 2023). While these effects will need to be accounted for
in future work, the tests presented here have demonstrated that the OC stream
contains a wealth of information that we have yet to understand.

Figure 24: Closest approach to the LMC. Top panel shows the stream on the sky,
coloured by the closest approach distance of each stream particle to the LMC.
Bottom panel shows the time of the closest approach. The closest approach
occurs at a distance of 5.8 kpc at a time of 366 Myr ago, at $\phi_{1}=-82.4$
deg. For the data, we fit in this work (i.e. $\phi_{1}\geq-70$ deg), the
closest approach is in the bin at $\phi_{1}=-70$ deg at a distance of 19.1 kpc
with a lookback time of 254 Myr. Figure 25: Uncertainty in the LMC’s orbit
over the past 500 Myr. The black curves show the LMC’s past orbit based on the
posterior samples in this work in Galactocentric coordinates. The red lines
show the LMC’s past orbit in the potential from McMillan (2017). The black and
red points show the LMC’s location at a lookback time of 254 Myr, i.e. the
point of closest approach to the OC stream. This shows that fits to the OC
stream can be used to meaningfully constrain the past location of the LMC.

###  6.6 Evolution in energy and angular momentum

Figure 26: The comparison of the energy vs angular momentum distribution in
the simulations vs data. The coloured pixel show the distribution of particles
from the best-fit model. The large circles show energy and angular momentum
measurements along the stream from this paper as presented in Figure 11. Each
point is coloured by $\phi_{1}$ similarly to Figure 11.

In order to further understand the LMC’s impact on the stream, as well as
understand the peculiar “V” shape seen in the observed energy versus angular
momentum (see Fig. 11), we now investigate the distribution of energy and
angular momentum in the model. In Figure 26, we show the observed and
simulated streams in the energy versus angular momentum plane. Note that in
order to aid this comparison, we use the same potential as we used in Section
4.3 (i.e., the best-fit model from McMillan, 2017) and not the time-varying
potential we have fit in this work. We see that the model reproduces the “V”
shape. This is not very surprising given the overall good agreement in
reproducing the phase-space track of the stream.

In order to understand how peculiar this distribution is, we first consider a
stream disrupting in a static potential. For such a stream, we would expect
the leading arm ($\phi_{1}>0$ deg) to have lower energies and lower absolute
angular momentum, and the trailing arm ($\phi_{1}<0$ deg) to have higher
energies and higher absolute angular momentum (see Gibbons et al., 2014).
Thus, we would broadly expect the energy and angular momentum distribution to
go from the upper left to the bottom right of Figure 26, with the smallest
$\phi_{1}$ in the upper left and the largest $\phi_{1}$ in the bottom right,
and $\phi_{1}$ smoothly varying in-between. Instead, the leading arm has more
absolute angular momentum and the same energy as the trailing arm, and the
trailing arm has the wrong correlation along the stream (i.e. $E$ increases as
${\rm L}_{z}$ increases). Furthermore, while the leading arm seemingly agrees
with the static model from Gibbons et al. (2014), $\phi_{1}$ should be
increasing as ${\rm L}_{z}$ increases: instead, it does the opposite.

With our best-fit model, we can understand this remarkable transformation.
From looking at different snapshots (see Appendix APPENDIX E:), we see that
the before the LMC’s infall, the stream was arranged as expected from the
static disruption model in Gibbons et al. (2014). During the LMC’s close
passage, the trailing arm is essentially flipped over in $E-L_{z}$ space,
while the leading arm has a substantial increase in $L_{z}$ and a decrease in
energy. This substantial re-arrangement of the stream in $E-L_{z}$ is similar
to the expected effect of the Sagittarius dwarf on stellar streams (Dillamore
et al., 2022).

##  7 Summary and Conclusions

In this paper, for the first time, using the combinations of Gaia, ${S}^{5}$,
and other spectroscopic survey data, we mapped the Orphan-Chenab stream over
200 degrees on the sky. Not only were we able to identify over 300
spectroscopic members associated with the OC stream, we also derived the full
6-D phase space track of the stream. We also measured the stream density along
the $\sim 120$ kpc of the extent of the structure. Here, we summarise the key
findings from this rich dataset:

  * •

We observe that while the apparent number density of stars along the stream
increases rapidly toward the stream edges, when corrected for the accordion
effect (slower speed at apocenter and faster speed at pericenter), the stream
density (or stream flow of stars) is mostly constant along the stream.

  * •

We do not detect a significant metallicity gradient in either the
Northern/Southern arms, or across the whole stream, preventing us from using
this to locate the current or previous position of the stream progenitor.

  * •

The mean metallicity of the stream based on ${S}^{5}$ measurements is $-1.9$
with a spread of 0.3 dex. Together with an estimate of the total stream
luminosity of $M_{V}=-8.2$, this makes the stream similar in properties to
classical MW satellites like Draco.

  * •

While the stream radial velocity dispersion is constant at around $\sim$ 5
$\mathrm{\,km\,s^{-1}}$ throughout the stream, we notice a significant change
in the intrinsic spread in proper motions, detecting a $\sim$ 10
$\mathrm{\,km\,s^{-1}}$ dispersion in velocities along the stream near its
pericenter.

  * •

The stream width changes significantly from 200 pc for the nearby part of the
stream to nearly 1 kpc for the most distant parts. This is accompanied by a
significant distance spread along the line of sight in the North, which has
been noticed previously by Sesar et al. (2013).

  * •

We detect large changes in the energy and the $z$-component of angular
momentum across the stream caused by the interaction with the LMC. This
spreading of the stream in the space of conserved quantities is expected to
significantly hamper our ability to use clustering in this space to identify
old accretion events.

  * •

The stars in the Southern part of the stream move at a significant angle with
respect to the stream track, both on the sky and along the line of sight.
Their velocity is clearly deflected in the direction of the LMC. This
deflection is the cause of the observed bending of the stream.

  * •

The model of the stream disruption in the presence of the LMC can successfully
reproduce the full 6-D stream track when the Galactic DM halo model is
represented by an oblate NFW-like potential. A prolate shape for the DM halo
provides a somewhat worse fit to the data, while a spherical shape is
inconsistent with the observations.

  * •

The inference of the MW dark matter halo shape has a degeneracy between
prolate and oblate shapes that can provide similar force fields in the OC
stream plane, and similar mass profiles in the radial distance range.

  * •

The modelling of the stream allows us to constrain the mass profile of both
the MW and the LMC, with the most precise estimates of enclosed masses for the
Milky Way and LMC being $M_{\mathrm{MW}}(<32.4\,\mathrm{kpc})=(2.85\pm
0.10)\times 10^{11}\mathrm{\,M_{\odot}}$,
$M_{\mathrm{LMC}}(<32.8\,\mathrm{kpc})=(7.02\pm 0.9)\times
10^{10}\,\mathrm{\,M_{\odot}}$, respectively. The best fit virial mass for the
MW is $M_{\mathrm{vir}}=7.74^{+4.14}_{-1.48}\times
10^{11}\mathrm{\,M_{\odot}}$ and the total mass for the LMC is
$M_{\mathrm{LMC}}=1.29^{+0.28}_{-0.23}\times 10^{11}\,\mathrm{\,M_{\odot}}$.

  * •

With the closest approach distance of the OC stream to the LMC of 21.3 kpc, we
can detect that the DM halo of the LMC extends out to at least 53 kpc.

  * •

The rewinding of MW satellites in the best-fit potential together with the OC
stream reveals that Grus II was likely accreted together with the OC
progenitor, as it has multiple close and low velocity encounters with the OC
stream.

In this paper, we have demonstrated the power of combining Gaia data with
dedicated spectroscopic follow-up by ${S}^{5}$ of this very long stellar
stream perturbed by the LMC. This has allowed us to put strong constraints on
the LMC and MW potential. It also shows that we can directly constrain the
past trajectory of the LMC and the amount of dynamical friction it
experienced, and potentially probe alternative gravity models that predict a
different past orbit of the LMC (e.g. MOND: Wu et al., 2008; Schee et al.,
2013).

The future of this analysis is manyfold. First, even the model presented in
this paper does not fully reproduce all of the observable quantities
perfectly, requiring a more complex model – such as a deforming and/or
triaxial MW DM halo – which our data should be able to constrain. Indeed,
Lilleengen et al. (2023) have already shown that the effect of the deforming
Milky Way on the OC stream is significant. Second, while in this work we used
an OC progenitor of fixed mass that was conveniently hidden in the part of the
stream located behind the Galactic disc, the properties of the progenitor –
such as its stellar and dark matter mass, and their extent – will likely have
an effect on observables, and can be constrained with our high-quality data.
With a detailed stream map and a stellar density measurement, we should also
be able to constrain the details of the mass loss history of the progenitor.

While the analysis presented here is potentially applicable to other streams
in the MW (e.g. Mateu, 2022), not many of them have detailed tracks (Patrick
et al., 2022), especially in 6-D (Li et al., 2022). This may however change
with the arrival of upcoming spectroscopic surveys like WEAVE (Jin et al.,
2022), DESI (Cooper et al., 2022) and 4MOST (de Jong et al., 2019). However,
the majority of streams are not as long as the Sgr or OC stream: thus,
limiting the amount of information we will be able to extract from each one.

In summary, despite the detailed analysis in this work, the OC stream still
contains many mysteries. Where is/was its progenitor? What were the properties
of that progenitor? What can it tell us about the Milky Way’s dark matter halo
and how it deforms in the presence of the LMC? We hope to explore these in
future work.

## Acknowledgements

SEK, DE, TSL, GFL, DZ & SM acknowledge funding through ARC DP210100855. TSL
acknowledges financial support from Natural Sciences and Engineering Research
Council of Canada (NSERC) through grant RGPIN-2022-04794. This research has
also been supported in part by the Australian Research Council Centre of
Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D), through
project number CE170100013. The paper includes data obtained with the Anglo-
Australian Telescope in Australia. We acknowledge the traditional owners of
the land on which the AAT stands, the Gamilaraay people, and pay our respects
to elders past, present and emerging.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

This project used public archival data from the Dark Energy Survey (DES).
Funding for the DES Projects has been provided by the U.S. Department of
Energy, the U.S. National Science Foundation, the Ministry of Science and
Education of Spain, the Science and Technology Facilities Council of the
United Kingdom, the Higher Education Funding Council for England, the National
Center for Supercomputing Applications at the University of Illinois at
Urbana-Champaign, the Kavli Institute of Cosmological Physics at the
University of Chicago, the Center for Cosmology and Astro-Particle Physics at
the Ohio State University, the Mitchell Institute for Fundamental Physics and
Astronomy at Texas A&M University, Financiadora de Estudos e Projetos,
Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro,
Conselho Nacional de Desenvolvimento Científico e Tecnológico and the
Ministério da Ciência, Tecnologia e Inovação, the Deutsche
Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy
Survey. The Collaborating Institutions are Argonne National Laboratory, the
University of California at Santa Cruz, the University of Cambridge, Centro de
Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the
University of Chicago, University College London, the DES-Brazil Consortium,
the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH)
Zürich, Fermi National Accelerator Laboratory, the University of Illinois at
Urbana-Champaign, the Institut de Ciències de l’Espai (IEEC/CSIC), the
Institut de Física d’Altes Energies, Lawrence Berkeley National Laboratory,
the Ludwig-Maximilians Universität München and the associated Excellence
Cluster Universe, the University of Michigan, the National Optical Astronomy
Observatory, the University of Nottingham, The Ohio State University, the
OzDES Membership Consortium, the University of Pennsylvania, the University of
Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the
University of Sussex, and Texas A&M University. Based in part on observations
at Cerro Tololo Inter-American Observatory, National Optical Astronomy
Observatory, which is operated by the Association of Universities for Research
in Astronomy (AURA) under a cooperative agreement with the National Science
Foundation.

The Legacy Surveys consist of three individual and complementary projects: the
Dark Energy Camera Legacy Survey (DECaLS; NOAO Proposal ID # 2014B-0404; PIs:
David Schlegel and Arjun Dey), the Beijing-Arizona Sky Survey (BASS; NOAO
Proposal ID # 2015A-0801; PIs: Zhou Xu and Xiaohui Fan), and the Mayall z-band
Legacy Survey (MzLS; NOAO Proposal ID # 2016A-0453; PI: Arjun Dey). DECaLS,
BASS and MzLS together include data obtained, respectively, at the Blanco
telescope, Cerro Tololo Inter-American Observatory, National Optical Astronomy
Observatory (NOAO); the Bok telescope, Steward Observatory, University of
Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOAO. The
Legacy Surveys project is honored to be permitted to conduct astronomical
research on Iolkam Du’ag (Kitt Peak), a mountain with particular significance
to the Tohono O’odham Nation.

NOAO is operated by the Association of Universities for Research in Astronomy
(AURA) under a cooperative agreement with the National Science Foundation.

The Legacy Surveys imaging of the DESI footprint is supported by the Director,
Office of Science, Office of High Energy Physics of the U.S. Department of
Energy under Contract No. DE-AC02-05CH1123, by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract; and by the U.S. National Science Foundation, Division of
Astronomical Sciences under Contract No. AST-0950945 to NOAO.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS-IV acknowledges support and resources from
the Center for High Performance Computing at the University of Utah. The SDSS
website is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic
Telescope LAMOST) is a National Major Scientific Project built by the Chinese
Academy of Sciences. Funding for the project has been provided by the National
Development and Reform Commission. LAMOST is operated and managed by the
National Astronomical Observatories, Chinese Academy of Sciences.

This paper made use of the Whole Sky Database (wsdb) created by Sergey Koposov
and maintained at the Institute of Astronomy, Cambridge by Sergey Koposov with
financial support from the Science & Technology Facilities Council (STFC) and
the European Research Council (ERC).

This paper used astropy (Astropy Collaboration et al., 2022), agama (Vasiliev,
2019), ipython (Pérez & Granger, 2007), imf https://github.com/keflavich/imf,
matplotlib (Hunter, 2007), minimint (Koposov, 2021), numpy (Harris et al.,
2020), scipy (Virtanen et al., 2020), sqlutilpy (Koposov, 2022b).

This paper made use of the Q3C package (Koposov & Bartunov, 2006).

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## Data Availability

This paper relies on publicly available data from Gaia, SDSS, APOGEE, LAMOST
surveys. We also use the data from S5 survey. The latest data is not yet
publicly available, but the previous release is available in Zenodo (Li & S5
Collaboration, 2021). The spectroscopic sample constructed in Section 3.2,
Stan models, measured parameters, the best fit OC stream model and posterior
samples for the stream models are provided either in the appendix or through
Zenodo https://zenodo.org/record/7222654.

## References

  * Abolfathi et al. (2018) Abolfathi B., et al., 2018, ApJS, 235, 42
  * Ahumada et al. (2020) Ahumada R., et al., 2020, ApJS, 249, 3
  * Amorisco et al. (2016) Amorisco N. C., Gómez F. A., Vegetti S., White S. D. M., 2016, MNRAS, 463, L17
  * Astropy Collaboration et al. (2022) Astropy Collaboration et al., 2022, ApJ, 935, 167
  * Baumgardt & Vasiliev (2021) Baumgardt H., Vasiliev E., 2021, MNRAS, 505, 5957
  * Belokurov et al. (2006) Belokurov V., et al., 2006, ApJ, 642, L137
  * Belokurov et al. (2007) Belokurov V., et al., 2007, ApJ, 658, 337
  * Belokurov et al. (2019) Belokurov V., Deason A. J., Erkal D., Koposov S. E., Carballo-Bello J. A., Smith M. C., Jethwa P., Navarrete C., 2019, MNRAS, 488, L47
  * Bennett & Bovy (2019) Bennett M., Bovy J., 2019, MNRAS, 482, 1417
  * Bernard et al. (2014) Bernard E. J., et al., 2014, MNRAS, 443, L84
  * Besla et al. (2007) Besla G., Kallivayalil N., Hernquist L., Robertson B., Cox T. J., van der Marel R. P., Alcock C., 2007, ApJ, 668, 949
  * Binney (2008) Binney J., 2008, MNRAS, 386, L47
  * Bonaca & Hogg (2018) Bonaca A., Hogg D. W., 2018, ApJ, 867, 101
  * Bonaca et al. (2021) Bonaca A., et al., 2021, ApJ, 909, L26
  * Boubert & Everall (2020) Boubert D., Everall A., 2020, MNRAS, 497, 4246
  * Bovy (2014) Bovy J., 2014, ApJ, 795, 95
  * Bovy (2015) Bovy J., 2015, ApJS, 216, 29
  * Bovy et al. (2016) Bovy J., Bahmanyar A., Fritz T. K., Kallivayalil N., 2016, ApJ, 833, 31
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Bullock & Johnston (2005) Bullock J. S., Johnston K. V., 2005, ApJ, 635, 931
  * Carlberg (2020) Carlberg R. G., 2020, ApJ, 889, 107
  * Carlberg & Grillmair (2013) Carlberg R. G., Grillmair C. J., 2013, ApJ, 768, 171
  * Carpenter et al. (2017) Carpenter B., et al., 2017, Journal of Statistical Software, 76, 1–32
  * Chabrier (2005) Chabrier G., 2005, in Corbelli E., Palla F., Zinnecker H., eds, Astrophysics and Space Science Library Vol. 327, The Initial Mass Function 50 Years Later. p. 41 (arXiv:astro-ph/0409465), doi:10.1007/978-1-4020-3407-7_5
  * Choi et al. (2016) Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., Johnson B. D., 2016, ApJ, 823, 102
  * Clementini et al. (2019) Clementini G., et al., 2019, A&A, 622, A60
  * Conroy et al. (2021) Conroy C., Naidu R. P., Garavito-Camargo N., Besla G., Zaritsky D., Bonaca A., Johnson B. D., 2021, Nature, 592, 534
  * Cooper et al. (2022) Cooper A. P., et al., 2022, arXiv e-prints, p. arXiv:2208.08514
  * Correa Magnus & Vasiliev (2022) Correa Magnus L., Vasiliev E., 2022, MNRAS, 511, 2610
  * Cui et al. (2012) Cui X.-Q., et al., 2012, Research in Astronomy and Astrophysics, 12, 1197
  * Cullinane et al. (2020) Cullinane L. R., et al., 2020, MNRAS, 497, 3055
  * de Boer et al. (2020) de Boer T. J. L., Erkal D., Gieles M., 2020, MNRAS, 494, 5315
  * de Jong et al. (2019) de Jong R. S., et al., 2019, The Messenger, 175, 3
  * Deason et al. (2021) Deason A. J., et al., 2021, MNRAS, 501, 5964
  * Dehnen & Binney (1998) Dehnen W., Binney J., 1998, MNRAS, 294, 429
  * Dillamore et al. (2022) Dillamore A. M., Belokurov V., Evans N. W., Price-Whelan A. M., 2022, MNRAS, 516, 1685
  * Dotter (2016) Dotter A., 2016, ApJS, 222, 8
  * Drimmel & Poggio (2018) Drimmel R., Poggio E., 2018, Research Notes of the American Astronomical Society, 2, 210
  * Erkal et al. (2016a) Erkal D., Sanders J. L., Belokurov V., 2016a, MNRAS, 461, 1590
  * Erkal et al. (2016b) Erkal D., Belokurov V., Bovy J., Sanders J. L., 2016b, MNRAS, 463, 102
  * Erkal et al. (2017) Erkal D., Koposov S. E., Belokurov V., 2017, MNRAS, 470, 60
  * Erkal et al. (2019) Erkal D., et al., 2019, MNRAS, 487, 2685
  * Erkal et al. (2021) Erkal D., et al., 2021, MNRAS, 506, 2677
  * Fardal et al. (2015) Fardal M. A., Huang S., Weinberg M. D., 2015, MNRAS, 452, 301
  * Fellhauer et al. (2007) Fellhauer M., et al., 2007, MNRAS, 375, 1171
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Gaia Collaboration et al. (2018) Gaia Collaboration et al., 2018, A&A, 616, A12
  * Gaia Collaboration et al. (2021) Gaia Collaboration et al., 2021, A&A, 649, A1
  * Garavito-Camargo et al. (2019) Garavito-Camargo N., Besla G., Laporte C. F. P., Johnston K. V., Gómez F. A., Watkins L. L., 2019, ApJ, 884, 51