# GD-1 Stellar Stream and Cocoon in the DESI Early Data Release

M. Valluri Department of Astronomy, University of Michigan, Ann Arbor, MI
48109, USA<EMAIL_ADDRESS>P. Fagrelius  NSF NOIRLab, 950 N. Cherry
Ave., Tucson, AZ 85719, USA  S. E. Koposov Institute for Astronomy, University
of Edinburgh, Royal Observatory, Blackford Hill, Edinburgh EH9 3HJ, UK
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3
0HA, UK  T. S. Li Department of Astronomy & Astrophysics, University of
Toronto, Toronto, ON M5S 3H4, Canada  Oleg Y. Gnedin Department of Astronomy,
University of Michigan, Ann Arbor, MI 48109, USA  E. F. Bell Department of
Astronomy, University of Michigan, Ann Arbor, MI 48109, USA  R. G. Carlberg
Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S
3H4, Canada  A. P. Cooper Institute of Astronomy and Department of Physics,
National Tsing Hua University, 101 Kuang-Fu Rd. Sec. 2, Hsinchu 30013, Taiwan
J. Aguilar  Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
CA 94720, USA  S. Ahlen Physics Dept., Boston University, 590 Commonwealth
Avenue, Boston, MA 02215, USA  C. Allende Prieto Departamento de Astrofísica,
Universidad de La Laguna (ULL), E-38206, La Laguna, Tenerife, Spain  Instituto
de Astrofísica de Canarias, C/ Vía Láctea, s/n, E-38205 La Laguna, Tenerife,
Spain  V. Belokurov  Institute of Astronomy, Madingley Rd, Cambridge CB3 0HA,
UK  L. Beraldo e Silva Steward Observatory, University of Arizona, 933 N,
Cherry Ave, Tucson, AZ 85721, USA  D. Brooks  Department of Physics &
Astronomy, University College London, Gower Street, London, WC1E 6BT, UK  A.
Byström Institute for Astronomy, University of Edinburgh, Royal Observatory,
Blackford Hill, Edinburgh EH9 3HJ, UK  T. Claybaugh  Lawrence Berkeley
National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  K. Dawson
Department of Physics and Astronomy, The University of Utah, 115 South 1400
East, Salt Lake City, UT 84112, USA  A. Dey NSF NOIRLab, 950 N. Cherry Ave.,
Tucson, AZ 85719, USA  P. Doel  Department of Physics & Astronomy, University
College London, Gower Street, London, WC1E 6BT, UK  J. E. Forero-Romero
Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio
Ip, CP 111711, Bogotá, Colombia  E. Gaztañaga  Institut d’Estudis Espacials de
Catalunya (IEEC), 08034 Barcelona, Spain  Institute of Cosmology and
Gravitation, University of Portsmouth, Dennis Sciama Building, Portsmouth, PO1
3FX, UK  Institute of Space Sciences, ICE-CSIC, Campus UAB, Carrer de Can
Magrans s/n, 08913 Bellaterra, Barcelona, Spain  S. Gontcho A Gontcho Lawrence
Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  Han,
J Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  K. Honscheid  Center for Cosmology and AstroParticle
Physics, The Ohio State University, 191 West Woodruff Avenue, Columbus, OH
43210, USA  Department of Physics, The Ohio State University, 191 West
Woodruff Avenue, Columbus, OH 43210, USA  The Ohio State University, Columbus,
43210 OH, USA  T. Kisner Lawrence Berkeley National Laboratory, 1 Cyclotron
Road, Berkeley, CA 94720, USA  A. Kremin Lawrence Berkeley National
Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  A. Lambert  Lawrence
Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  M.
Landriau  Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
CA 94720, USA  L. Le Guillou Sorbonne Université, CNRS/IN2P3, Laboratoire de
Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France  M.
E. Levi Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA
94720, USA  A. de la Macorra Instituto de Física, Universidad Nacional
Autónoma de México, Cd. de México C.P. 04510, México  M. Manera Departament de
Física, Serra Húnter, Universitat Autònoma de Barcelona, 08193 Bellaterra
(Barcelona), Spain  Institut de Física dâ€™Altes Energies (IFAE), The
Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra
Barcelona, Spain  P. Martini Center for Cosmology and AstroParticle Physics,
The Ohio State University, 191 West Woodruff Avenue, Columbus, OH 43210, USA
Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory,
140 W 18th Avenue, Columbus, OH 43210, USA  The Ohio State University,
Columbus, 43210 OH, USA  G. E. Medina Department of Astronomy & Astrophysics,
University of Toronto, Toronto, ON M5S 3H4, Canada  A. Meisner NSF NOIRLab,
950 N. Cherry Ave., Tucson, AZ 85719, USA  R. Miquel  Institució Catalana de
Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010 Barcelona,
Spain  Institut de Física d’Altes Energies (IFAE), The Barcelona Institute of
Science and Technology, Campus UAB, 08193 Bellaterra Barcelona, Spain  J.
Moustakas Department of Physics and Astronomy, Siena College, 515 Loudon Road,
Loudonville, NY 12211, USA  A. D. Myers  Department of Physics & Astronomy,
University of Wyoming, 1000 E. University, Dept. 3905, Laramie, WY 82071, USA
J. Najita NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA  C. Poppett
Space Sciences Laboratory, University of California, Berkeley, 7 Gauss Way,
Berkeley, CA 94720, USA  University of California, Berkeley, 110 Sproul Hall
#5800 Berkeley, CA 94720, USA  F. Prada Instituto de Astrofísica de Andalucía
(CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain  M. Rezaie
Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan,
KS 66506, USA  G. Rossi  Department of Physics and Astronomy, Sejong
University, Seoul, 143-747, Korea  A. H. Riley Institute for Computational
Cosmology, Department of Physics, Durham University, South Road, Durham DH1
3LE, UK  E. Sanchez CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain  D.
Schlegel  Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley,
CA 94720, USA  M. Schubnell  Department of Physics, University of Michigan,
Ann Arbor, MI 48109, USA  University of Michigan, Ann Arbor, MI 48109, USA  D.
Sprayberry  NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA  G. Tarlé
University of Michigan, Ann Arbor, MI 48109, USA  G. Thomas Instituto de
Astrofísica de Canarias, C/ Vía Láctea, s/n, E-38205 La Laguna, Tenerife,
Spain Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206, La
Laguna, Tenerife, Spain  B. A. Weaver  NSF NOIRLab, 950 N. Cherry Ave.,
Tucson, AZ 85719, USA  R. H. Wechsler Kavli Institute for Particle
Astrophysics and Cosmology, Stanford University, Menlo Park, CA 94305, USA
Physics Department, Stanford University, Stanford, CA 93405, USA  SLAC
National Accelerator Laboratory, Menlo Park, CA 94305, USA  R. Zhou Lawrence
Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA  H.
Zou National Astronomical Observatories, Chinese Academy of Sciences, A20
Datun Rd., Chaoyang District, Beijing, 100012, P.R. China

###### Abstract

We present $\sim 126$ new spectroscopically identified members of the GD-1
tidal stream obtained with the 5000-fiber Dark Energy Spectroscopic Instrument
(DESI). We confirm the existence of a “cocoon” which is broad (FWHM$\sim
2.932^{\circ}\sim 460$pc) and kinematically hot (velocity dispersion,
$\sigma\sim 5-8$$\,\mathrm{km\ s}^{-1}$) component that surrounds a narrower
(FWHM$\sim 0.353^{\circ}\sim 55$pc) and colder ($\sigma\sim
2.2-2.6$$\,\mathrm{km\ s}^{-1}$) thin stream component (based on a median per
star velocity precision of 2.7 $\,\mathrm{km\ s}^{-1}$). The cocoon extends
over at least a $\sim 20^{\circ}$ segment of the stream observed by DESI. The
thin and cocoon components have similar mean values of [Fe/H]: $-2.54\pm 0.04$
dex and $-2.45\pm 0.06$ dex suggestive of a common origin. The data are
consistent with the following scenarios for the origin of the cocoon. The
progenitor of the GD-1 stream was an accreted globular cluster (GC) and: (a)
the cocoon was produced by pre-accretion tidal stripping of the GC while it
was still inside its parent dwarf galaxy; (b) the cocoon is debris from the
parent dwarf galaxy; (c) an initially thin GC tidal stream was heated by
impacts from dark subhalos in the Milky Way; (d) an initially thin GC stream
was heated by a massive Sagittarius dwarf galaxy; or a combination of some
these. In the first two cases the velocity dispersion and mean metallicity are
consistent with the parent dwarf galaxy having a halo mass of $\sim
10^{9}$$\,M_{\odot}$. Future DESI spectroscopy and detailed modeling may
enable us to distinguish between these possible origins.

Unified Astronomy Thesaurus concepts: Atomic spectroscopy (2099), Spectroscopy
(1558), Dark matter (353), Stellar streams (2166), Milky Way dynamics (1051);
Milky Way dark matter halo (1049)

††journal: ApJ††software: astropy (Astropy Collaboration et al., 2018, 2022),
numpy (Harris et al., 2020), scipy (Virtanen et al., 2020), scikit-learn
(Pedregosa et al., 2011), matplotlib (Hunter, 2007), emcee (Foreman-Mackey et
al., 2013), gala (Price-Whelan, 2017; Price-Whelan et al., 2022).
††facilities: KPNO: Mayall (DESI), Gaia, Blanco (DECam).

##  1 Introduction

Stellar streams are a stunning example of the ongoing hierarchical assembly of
the Milky Way via the gravitational interaction between our Galaxy and its
satellites. They are the result of the disruptive tidal forces of the Galaxy
on its companions such as globular clusters and dwarf galaxies. If the
progenitor of a tidal stream is of low enough mass, the tidal debris that is
liberated follows a path that is close to the orbit of its progenitor.

Long tidal streams that extend over tens or even hundreds of degrees on the
sky are of particular interest since they can be used to probe the global mass
density profile, total mass and 3D shape of the Milky Way’s dark matter halo
(e.g., Helmi, 2004; Johnston et al., 2005; Law et al., 2009; Law & Majewski,
2010; Koposov et al., 2010; Malhan & Ibata, 2019; Vasiliev et al., 2021). A
stellar stream is the result of tidal fields removing stars from the outskirts
of a progenitor at a velocity dispersion close to the cluster velocity
dispersion, with the tide at the stream pericenter creating an outward ordered
angular momentum distribution that leads to orbital shearing and a wider
orbital spread at apocenter. Since globular clusters (GCs) have small internal
velocity dispersions ($\lesssim 2$$\,\mathrm{km\ s}^{-1}$), unperturbed tidal
streams produced by GCs are both thin and dynamically cold. Therefore velocity
dispersion variations, density perturbations and gaps along a GC stream might
hold clues to past fly-by interactions and impact from dark matter subhalos
which can produce gaps and heating which results in broadening the streams
over time. If the streams do not pass through the Galactic disk they may be
especially good antennae for detecting such dark subhalos. Consequently
observations of both the spatial and velocity substructure in tidal streams
have elicited great interest as potential probes of dark matter substructure
(e.g., Johnston et al., 2002; Ibata et al., 2002; Siegal-Gaskins & Valluri,
2008; Carlberg, 2012, 2013; Erkal & Belokurov, 2015; Bonaca et al., 2019).

The GD-1 stellar stream, first identified by Grillmair & Dionatos (2006),
likely originated as a GC (Koposov et al., 2010). It is located in the
Galactic halo, and being both thin and very long (spanning at least 100∘), is
a good candidate for exploring both the global distribution of dark matter in
the inner Milky Way (e.g., Koposov et al., 2010; Bowden et al., 2015; Malhan &
Ibata, 2019) and for detecting dark subhalos from the analysis of gaps and
velocity dispersion variations. The progenitor remains undetected, possibly
having been completely disrupted (de Boer et al., 2018; Malhan et al., 2018;
Price-Whelan & Bonaca, 2018). GD-1 is of great current interest since Gaia DR2
data led to the detection of density variations, gaps and off-stream
structures e.g., a ‘spur’ a ‘blob’ and ‘wiggles’ (Price-Whelan & Bonaca, 2018;
de Boer et al., 2020) which have been claimed as evidence of possible impacts
by dark matter subhalos (Bonaca et al., 2019). Since the stream does not pass
through the inner disk region (it has a perigalactic (apogalactic) distance of
$\sim 14\leavevmode\nobreak\ (26)$ kpc, and inclination of 39∘ to the disk
plane, Koposov et al. (2010)) it is unlikely to have been heated by encounters
with giant molecular clouds in the inner disk. It has also been shown that
GD-1’s orbit has a very low probability of impact with other (intact) GCs
(Doke & Hattori, 2022).

An alternative explanation for the gaps and density variations in the GD-1
stream has been proposed: regularly spaced gaps have been predicted to arise
naturally due to epicyclic motions of stars along the stream (Küpper et al.,
2012, 2015). From an analysis of the spectrum of stream gaps and after
accounting for artifacts in the Gaia DR2 scanning law and projection effects,
Ibata et al. (2020) conclude that the density profile of the stream shows
periodic gaps separated by $\sim 2.64\pm 0.18$kpc and that such a regularly
spaced gap distribution can be produced via epicyclic motions in a smooth
Galactic potential. However since the regularly spaced gaps arising from
internal dynamics are not predicted to be associated with increased velocity
dispersions, obtaining kinematic data is useful for distinguishing this
process from the effects of perturbations.

Malhan et al. (2019b) used Gaia DR2 data to show that the GD-1 stellar stream
possesses a secondary diffuse and extended stellar component ($\sim 100$ pc
wide) that surrounds the previously identified thinner component of the
stream. They named this broader component the “cocoon” and showed that it was
detected at $>5\sigma$ confidence level. Similar broad and extended
morphological and kinematic features have also been found in several other
streams (e.g., Jhelum, Phlegethon, Malhan et al., 2020, 2022a; Ibata et al.,
2023; Awad et al., 2024). Simulated GC streams that form in dwarf galaxies and
are then accreted into Milky Way like halos show complex morphological
features (like gaps, “spurs”, “blobs”, “cocoons”). Such accreted GC streams
have complex morphologies in both cosmological simulations (Carlberg, 2018;
Carlberg & Agler, 2023) and controlled simulations (Malhan et al., 2020).
These complex morphologies are produced via the tidal interactions between the
GCs and the central dark matter distributions in the GC’s parent dwarf galaxy
prior to their infall and subsequent accretion onto a massive host galaxy like
the Milky Way. Furthermore it has been shown that the phase space structure of
accreted GC streams could be used to obtain constraints on the dark matter
distribution in the parent dwarf galaxy (Malhan et al., 2020, 2022a), making
the phase space structure of the cocoon important to determine.

Recently (Dillamore et al., 2022a) showed that the GD-1 stream is extremely
sensitive to perturbations from the Sagittarius dwarf galaxy (Sgr) especially
if it was quite massive ($\gtrsim 4\times 10^{10}$$\,M_{\odot}$) since the
orbit of the GD-1 stream intersects with that of Sgr. They showed that a wide
range of outcomes is possible with various types of substructures produced in
and around the GD-1 stream (see their Fig 9). Interactions with the Sgr dwarf
change the energy distribution of the GD-1 particles making it substantially
broader; also the streams density profile is altered (affecting the surface
brightness distribution along the stream) and its length is modified. While
many of the induced perturbations are difficult to distinguish from the cocoon
formation scenario above, the change in length is particularly interesting as
this is a unique prediction of the heavy Sgr hypothesis.

Since impacts from dark matter sub-halos which cross a stream also cause
stream heating (in addition to producing shortlived gaps and off-stream
features like ‘spurs’), it is particularly important to measure the velocity
dispersions of robustly identified stream candidates. Thus the phase space
structure of the GD-1 stream has the potential to improve our understanding of
the nature of dark matter both within the Milky Way and in dwarf satellites
that were accreted by our Galaxy.

While Gaia photometry and proper-motions led to the detection of an extension
to the stream as well as these various complex morphological features, there
are fewer than 100 spectroscopically and chemically confirmed members from
previous surveys with high radial velocity uncertainty of $\sim
1-2$$\,\mathrm{km\ s}^{-1}$ (Bonaca et al., 2020; Gialluca et al., 2021; Ibata
et al., 2020). An additional $\sim 200$ spectroscopically confirmed members
are available from the SDSS and LAMOST surveys but have lower radial velocity
precision ($\sim 5-20$$\,\mathrm{km\ s}^{-1}$). In this work we present 126
new spectroscopically confirmed members of the GD-1 stream obtained with the
Dark Energy Spectroscopic Instrument (DESI) with $\sim 1-10$$\,\mathrm{km\
s}^{-1}$ velocity precision (median precision 2.7$\,\mathrm{km\ s}^{-1}$) and
$\sim 0.2$ dex metallicity precision. We also improve on 17 previous radial
velocity measurements.

Starting in December 2020 DESI carried out an approximately 6 month long
“Survey Validation” (SV) phase. The data from SV were used to validate
scientific requirements, specifically the target selection algorithms (DESI
Collaboration et al., 2022). During the SV period, while testing and
streamlining targeting algorithms and analysis pipelines, there was also an
opportunity to observe specific targets. In this period, DESI took several
observations in the GD-1 stream region. Since these data were obtained
primarily as part of the SV program, the target selection criteria and
exposure times were slightly different from the main survey. After the
completion of SV, there have continued to be observations in the GD-1 stream
area.

In this paper we introduce the spectroscopically identified stars associated
with the GD-1 tidal stream from the DESI SV observations, which was recently
released to the public (DESI Collaboration et al., 2023a; Koposov et al.,
2024). Additional data releases are planned at regular intervals. In section
2, we briefly describe the DESI Survey Validation data set in the region
overlapping GD-1. In section 3 we discuss how we identify GD-1 stream members
from the DESI targets using cuts in the stellar color-magnitude diagram, Gaia
parallax, Gaia proper-motions as well as radial velocity and metallicity
[Fe/H]. In section 4 we discuss the phase space distribution of stars in the
stream, focusing on the identification and kinematics of the broader region
referred to as the cocoon. Finally, in Section 5 we summarize our findings and
discuss the implications of these results. We also give a preview of what to
expect from future DESI data releases on the GD-1 stream region.

##  2 DESI observations of the GD-1 region

Figure 1: GD-1 stream in great-circle coordinates ($\phi_{1},\phi_{2}$) with
the stream along $\phi_{1}$. Grey scale shows 2D histogram of surface density
of Gaia DR3 stars with a broad color-magnitude, and a proper motion selection
as indicated in the legend. The GD-1 stream is clearly visible as an over-
density at roughly $\phi_{2}=0$. Magenta curves show the DESI tiles observed
during commissioning and Survey Validation.

The Dark Energy Spectroscopic Instrument (DESI) is a robotic, fiber-fed,
highly multiplexed spectrograph installed on the Mayall 4-meter telescope at
Kitt Peak National Observatory (DESI Collaboration et al., 2022). A central
feature of the instrument is the use of 5,000 robotic positioners that move
individual optical fibers to pre-identified targets over a $\sim 3^{\circ}$
field of view (DESI Collaboration et al., 2016; Silber et al., 2023; Miller et
al., 2023). The DESI survey is a five-year survey of about 14,000 square
degrees to obtain spectra for approximately 40 million galaxies and quasars
and over 10 million Milky Way stars (DESI Collaboration et al., 2016).

DESI was installed in 2019 and commissioning and survey validation extended
through May 2021 (DESI Collaboration et al., 2023b). Since May 2021, the
5-year DESI survey has been underway and has obtained an unprecedented number
of galaxy, quasar and stellar spectra. The primary scientific goal of the
survey is to use spectroscopic redshifts of galaxies and quasars to obtain the
most precise constraints on the expansion history of the universe ever
obtained using the baryon acoustic oscillation and other methods (e.g.,
redshift space distortions) (DESI Collaboration et al., 2016; Levi et al.,
2013).

Over the first 5-years of the DESI Milky Way Survey (Allende Prieto et al.,
2020; Cooper et al., 2023) it is expected that spectra of $\simeq 7$ million
unique Milky Way stars at Galactic latitudes $|b|>20^{\circ}$ will be obtained
during the main survey and millions more will be obtained during poor sky
conditions (the “back up program”). The DESI-Milky Way Survey (MWS) uses a
simple and inclusive target selection scheme focused on the thick disk and
stellar halo. All targets are selected from the Gaia catalog (originally DR2,
but subsequently DR3). Targets fall into three categories: MAIN-BLUE, MAIN-
RED, and MAIN-BROAD, within the magnitude range $16<r<19$. MAIN-BLUE
prioritizes all Gaia point sources in this magnitude range with blue optical
colors ($g-r<0.7$) in the DESI Legacy Imaging Surveys (Dey et al., 2019). The
MAIN-RED selection applies Gaia proper-motion and parallax criteria to sources
with redder colors $(g-r>0.7)$ to increase the probability of observing
distant halo giants. Sources with $g-r>0.7$ that do not meet the astrometric
criteria are targeted at lower priority in the MAIN-BROAD category (Cooper et
al., 2023). DESI-MWS is expected to achieve an average spectroscopic
completeness of $\sim 28$% over the various MAIN target classes.

The DESI-MWS is executed during bright-sky conditions (when high-redshift
galaxy observations are inefficient) and shares the focal plane with a low-
redshift Bright Galaxy Survey (BGS; Hahn et al., 2023), although BGS targets
have higher priority for fiber assignment than MWS stars.

The DESI spectroscopic data presented in this paper were taken during the on-
sky DESI SV from 2020/12/14 to 2021/5/14; the corresponding data (both spectra
and the catalogs used in this paper) are now publically available under DESI’s
Early Data Release (EDR) https://data.desi.lbl.gov/doc/releases/edr/vac/mws/.
During the SV, the DESI Collaboration carried out observations to evaluate the
targeting pipelines, instrument and observatory readiness in three separate
subprograms SV1, SV2 and SV3. The MWS target selection for SV was previously
described in Allende Prieto et al. (2020). Other DESI SV programs are
discussed in DESI Collaboration et al. (2023b). Each fiber configuration in a
given field is referred to as a “tile”. Some of these tiles had a larger
percentage of fibers dedicated to Milky Way stars than expected for the main
survey and in some cases went to fainter magnitudes (for the SV, stars in the
range $19<r<20$ were also included at lower priority.) As discussed in Section
7 of Cooper et al. (2023), several fields in the GD-1 region were observed
during SV1 but without any special prioritization of potential GD-1 targets.
In addition, 22 tiles were observed during SV2, also without any special
prioritization for GD-1 stars. In total 45 tiles in the GD-1 region were
observed during Survey Validation. For more details on the observing
conditions and various programs for SV, readers are directed to the MWS EDR
Value Added Catalogue paper (Koposov et al., 2024).

Figure 1 shows the GD-1 stream region in a coordinate system in which
$\phi_{1}$ is the angle on the sky along the stream and $\phi_{2}$ is the
angle perpendicular to $\phi_{1}$. The transformation from Gaia DR3 sky
positions ($\alpha$, $\delta$) and proper motions
($\mu_{\alpha}^{*},\mu_{\delta}$) to this natural GD-1 frame was carried out
using the coordinate transformation described in Koposov et al. (2010)111The
pole of the great circle frame fitting the GD-1 stream has (RA, Dec)=
(34.5987, 29.7331) in 2000 coordinates., which also yields proper motion in
[$\phi_{1},\phi_{2}$], hereafter referred to as $\mu_{\phi_{1}}$ and
$\mu_{\phi_{2}}$222Hereafter we drop the ∗ and simply use $\mu_{\alpha}$ and
$\mu_{\phi_{1}}$ to denote $\mu_{\alpha}\cos(\delta)$ and
$\mu_{\phi_{1}}\cos(\phi_{2})$ respectively. Most targets were also identified
in the DESI Legacy Survey (DECaLs) (Dey et al., 2019), from which the $g$ and
$r$-band magnitudes for these stars were measured. The Legacy Survey
magnitudes were de-reddened using dust maps (Schlafly et al., 2010) updated
for the Legacy Survey and available on the Legacy Survey website
https://www.legacysurvey.org. These data were cross-matched with the Gaia DR3
catalog for parallax and proper motion measurements (Gaia Collaboration,
2022). The Gaia DR3 proper motions were corrected for solar reflex motion
using Astropy. We used Astropy’s default parameters for
galactocentric_frame_defaults.set(‘v4.0’) (Reid & Brunthaler, 2004; GRAVITY
Collaboration et al., 2018; Drimmel & Poggio, 2018; Bennett & Bovy, 2019).
Distances to individual stars were estimated using the distance modulus in
Equation 1 below (also see Appendix). Gaia parallaxes were zero-point
corrected following the prescriptions in (Lindegren et al., 2021a).

The 2D histogram in Figure 1 shows the surface density of Gaia DR3 stars with
DECaLS $0.2\leq(g-r)<0.4$ and $18\leq r<20$, $-10$mas
yr${}^{-1}<\mu_{\phi_{1}}<-4$mas yr-1 and $|\mu_{\phi_{2}}|<-5$mas yr-1. These
selections are similar to and based on previous works (e.g., Price-Whelan &
Bonaca, 2018). The boundaries of the $45$ DESI tiles in this region that were
observed during SV are shown by magenta curves. The direction of motion of the
stream is from positive $\phi_{1}$ to negative $\phi_{1}$.

Figure 2: Gaia-DR3 stars in the GD-1 region defined in the frame of reference
of the stream [$|\phi_{2}|<3^{\circ}$, $-85^{\circ}<\phi_{1}<20^{\circ}$] with
a broad color-magnitude cut of $0.2\leq g-r<0.4$ and $18\leq r<20$. (a) 2D
histogram of $\mu_{\phi_{1}}$ vs. $\phi_{1}$ with a spline (red curve)
highlighting the GD-1 stream track. (b) 2D histogram of $\mu_{\phi_{2}}$ vs.
$\phi_{1}$ with $\mu_{\phi_{1}}$ within 2mas/yr from red stream track in
panel-(a). Cyan curve shows spline fit to stream. (c) scatter plot of
$\mu_{\phi_{2}}$ vs. $\mu_{\phi_{1}}$. Grey dots are all stars in panel-(b).
Black points are stars with separation from red and cyan splines with
$\delta\mu_{\phi_{1}}<2$mas/yr and $\delta\mu_{\phi_{2}}<2$mas/yr and parallax
and broad color-magnitude cuts as described in the text. (d) Selection in
panel-(c) shown as a 2-D histogram in $\phi_{1}-\phi_{2}$ (relative density in
arbitrary units) with a spline representing stream track shown in red. Figure
3: Color absolute magnitude diagram with the dereddened magnitudes from
DECaLS. For this plot stars were selected with $|\delta\phi_{2}|<0.2$deg from
the stream track and both $\mu_{\phi_{1}}$ and $\mu_{\phi_{2}}$ within
0.5mas/yr from the respective splines in Figure 2. This selection is only used
to define the color-magnitude selection. We use the distance modulus vs
$\phi_{1}$ relation in Eqn. 1 to correct for distance gradient along the
stream. Overlaid is an isochrone from the PARSEC Isochrone database with a
metallicity [M/H] = -2.0dex and age = 12Gyr (red curve). The light red curves
are at $\Delta|g-r|=0.1$ from the red curve.

The stellar spectra obtained by the DESI spectrograph are processed with the
RVS pipeline, based on RVSpec (Koposov, 2019), to measure each star’s line-of-
sight (radial) velocity and stellar atmospheric parameters. The spectra are
further processed by the SP pipeline based on FERRE (Allende Prieto et al.,
2006)333 https://github.com/callendeprieto/ferre and supporting python code
piferre https://github.com/callendeprieto/piferre, which measures [Fe/H] and
[$\alpha$/Fe], among other properties. The RVS pipeline measures the radial
velocity with a precision of $\leq 1-10$$\,\mathrm{km\ s}^{-1}$and [Fe/H] with
a precision of $\sim$0.15-0.2 dex. For further details of the pipelines and
their validation see Cooper et al. (2023) and Koposov et al. (2024).

##  3 GD-1 Stream Membership

###  3.1 Initial selection

Our search for GD-1 stream members begins by restricting to the region
$|\phi_{2}|<3^{\circ}$, $-85^{\circ}<\phi_{1}<20^{\circ}$ (galactic
coordinates: $210^{\circ}<l<90^{\circ}$, $20^{\circ}<b<70^{\circ}$). We also
apply the same broad cut in color-magnitude space and proper-motion selections
as used in Figure 1. These selections expand slightly on the limits set by
Price-Whelan & Bonaca (2018) and Malhan & Ibata (2019).

The distance to the GD-1 stream is known to vary with $\phi_{1}$ (Koposov et
al., 2010). We assume a distance modulus $D$ given by,

$\displaystyle D=18.82+[(\phi_{1}+50)/64]^{2}-4.45,$ (1)

which is similar to the function estimated by Koposov et al. (2010). This
relation is a simple quadratic fit to the subgiant branch stars in GD-1. More
details are provided in the appendix where we also compare this relation with
the distance relations proposed by several other authors (Price-Whelan &
Bonaca, 2018; de Boer et al., 2018; Li et al., 2018). The functions defined by
de Boer et al. (2018) and Li et al. (2018) differ only slightly from the
distance-$\phi_{1}$ relation in Equation 1 (see Figure 13 in Appendix A. We
use this equation in the rest of the paper. We impose a Gaia parallax cut such
that

$\displaystyle\delta\varpi$ $\displaystyle=|1/d(\phi_{1})-\varpi|,$ (2)
$\displaystyle\delta\varpi$ $\displaystyle<\varpi_{sys}+2*\varpi_{err},$ (3)

where $d(\phi_{1})$ is the distance to the stream at a given value $\phi_{1}$
computed assuming the distance modulus in Equation 1, $\varpi_{err}$ is the
uncertainty on the Gaia parallax for each star, and $\varpi_{sys}=0.05$mas is
the systematic error on the parallax (Lindegren et al., 2021b).

Figure 2 shows three projections of Gaia stars in the GD-1 stream region: the
stream is clearly offset from the background in $\mu_{\phi_{1}}$ vs.
$\phi_{1}$ (panel-a) and $\mu_{\phi_{2}}$ vs. $\phi_{1}$ (panel-b). We fit
splines to the stream tracks that emerge in panels a and b (shown by red and
cyan curves). We include all stars $\pm$ 2 mas/yr around each spline to obtain
the selection shown by black points in Figure 2 (c) which is similar to the
polygon defined in Price-Whelan & Bonaca (2018). Figure 2 (d) shows the stream
in $\phi_{1}-\phi_{2}$ coordinates over-plotted with a cubic spline through
the highest density regions of the stream shown in red. The values at the
knots as a function of $\phi_{1}$ for these three splines as well as for the
velocity spline in Figure 4 (a) (to be discussed in Section 3.3) are given in
Table 1.

Table 1: Knots for four stream-track splines
---
$\phi_{1}$ | $\mu_{\phi_{1}}$ | $\mu_{\phi_{2}}$ | ${\phi_{2}}$ | $V_{GSR}$
[∘] | [mas/yr] | [mas/yr] | [∘] | [km/s]
90.0 | -3.00 | -5.26 | -0.70 | 79.79
-70.0 | -1.50 | -7.08 | -0.87 | 59.15
-60.0 | -0.72 | -7.44 | -0.70 | 27.61
-46.0 | -0.10 | -7.43 | -0.38 | -21.66
-37.6 | 0.05 | -7.18 | -0.21 | -46.38
-29.2 | 0.22 | -6.78 | -0.09 | -66.70
-20.8 | 0.21 | -6.27 | 0.00 | -93.19
-12.4 | -0.18 | -5.66 | 0.08 | -113.29
-4.0 | -0.84 | -4.98 | 0.16 | -131.91
-15.0 | -0.02 | -5.85 | 0.06 | -107.81
0.0 | -1.20 | -4.64 | 0.20 | -144.31
20.0 | -2.50 | -2.90 | 0.50 | -257.76

###  3.2 Color-Magnitude Diagram selection

Given the above initial cuts as an input, we show that the data are well
identified in a color-magnitude diagram (CMD) as shown in Figure 3. The stars
are well matched to an isochrone of 12 Gyr with [M/H] $=-2.0$ from the PARSEC
isochrone database (Bressan et al., 2012; Tang et al., 2014; Marigo et al.,
2017; Pastorelli et al., 2019,
2020)444https://people.sissa.it/~sbressan/parsec.html. Other authors have used
similar isochrones selections (e.g., Ibata et al., 2020). We accepted stars
along the PARSEC isochrone if they lie within $\delta(g-r)\pm 0.1$ (shown by
faint red curves).

###  3.3 Radial Velocity Distribution

Figure 4: (a) Line-of-sight velocities (corrected for solar reflex motion and
transformed to GSR frame) for stars from DESI (black points), SDSS (Ahn et
al., 2012) or LAMOST (Yan et al., 2022) (blue diamonds), MMT-Hectochelle
(Bonaca et al., 2020) (orange triangles). All stars were selected with the
parallax, proper motion, and CMD cuts discussed in the text and include only
stars with $|\delta\phi_{2}|<3^{\circ}$. The stream is clearly visible and
highlighted by the red curve which is a spline fit (see text). All DESI radial
velocities and [Fe/H] measurements in this figure are given in Table 4 of
Appendix B. (b) Kernel density distributions of stars with velocities within
$\pm 30$$\,\mathrm{km\ s}^{-1}$ (i.e. between dashed red lines) from the red
stream track in left panel ($-85^{\circ}<\phi_{1}<20^{\circ}$). The kernel
band widths are 2.74$\,\mathrm{km\ s}^{-1}$ for DESI (black curve),
7.25$\,\mathrm{km\ s}^{-1}$ for SDSS+LAMOST (blue dot-dashed curve) and
1.5$\,\mathrm{km\ s}^{-1}$ for MMT (orange dashed curve).

With the stream defined by the cuts discussed in the previous subsections, we
now examine the radial velocities of the selected stars. Figure 4 (a) shows
the solar reflex-motion corrected line-of-sight velocities in the Galactic
Standard of Rest frame $V_{GSR}$ as a function of $\phi_{1}$ for stars in the
GD-1 region selected by parallax, proper-motion (as shown in Figure 2), CMD
cuts discussed in the text and include only stars with $|\delta\phi_{2}|\leq
3^{\circ}$. Stars observed by DESI from the EDR catalog (DESI Collaboration et
al., 2023a) are shown as black points. Stars observed by SDSS (Ahn et al.,
2012) (blue) DR9555https://www.sdss3.org/dr9/ and LAMOST DR8
666http://www.lamost.org/dr8/ are shown by blue symbols. The 43 stars
identified as GD-1 members using MMT-Hectochelle spectra (Bonaca et al., 2020)
are shown by orange symbols. The selection of (Bonaca et al., 2020) only
considered stars with velocities within 7$\,\mathrm{km\ s}^{-1}$ from the
stream track or the spur as GD-1 members, but we do not expand on this
selection since none of their stars belong to the cocoon region we discuss
later.

All DESI-EDR measurement in Figure 4 (a) are given in a machine readable FITS
file available at https://zenodo.org/records/11638330, the first 3 rows of the
table are shown in Table 4 in Appendix B.

We note that 18 stars from the SDSS and LAMOST datasets were also observed by
DESI. Where there were duplicates we selected the measurement with the smaller
velocity uncertainty. Except for one star DESI velocity uncertainties were
always smaller than uncertainties from SDSS and LAMOST. There were no overlaps
with GD-1 members in the MMT sample of Bonaca et al. (2020) primarily because
the DESI-EDR observations did not have much coverage of the spur area. Of the
remaining 153 SDSS+LAMOST stars 14 did not have [Fe/H] measurements and were
therefore excluded.

In summary DESI has added 126 new spectroscopically confirmed radial
velocities and improved RV measurements for 17 previously observed members. A
spline fit to the stream is shown in red (see below for details of how this
fit is obtained) along with dashed lines showing $\delta V_{GSR}=\pm
30$$\,\mathrm{km\ s}^{-1}$ from the stream track. Figure 4 (b) shows kernel
density estimates of $\delta V_{GSR}$ (from the stream track in the left
panel) of stars observed by DESI (black), SDSS+LAMOST (blue) and MMT (orange).
All stars within $-85^{\circ}<\phi_{1}<20^{\circ}$ are shown. A clear peak is
visible identifying the GD-1 stream. In making the kernel density estimate we
have used band-widths 2.74$\,\mathrm{km\ s}^{-1}$ for DESI (black curve),
7.25$\,\mathrm{km\ s}^{-1}$ for SDSS+LAMOST (blue dot-dashed curve) and
1.5$\,\mathrm{km\ s}^{-1}$ for MMT (orange dashed curve). The broader
histogram corresponding to stars from SDSS+LAMOST reflects their larger
measurement uncertainties, while the narrower MMT distribution reflects both
the small uncertainties and the more restrictive selection of $|\delta
V_{GSR}|\leq 7$$\,\mathrm{km\ s}^{-1}$ imposed by Bonaca et al. (2020).

In order to measure the velocity dispersion of the stream we need to properly
account for the uncertainties in the radial velocity measurements and account
for background contamination. We use a mixture modeling approach to determine
the velocity track and velocity dispersion of the stream. We use the approach
previously adopted for other streams in Koposov et al. (2019, 2023). Focusing
on the region along the stream of $-55^{\circ}\leq\phi_{1}\leq 0^{\circ}$,
where the density of the stream and number of DESI RV measurements is the
highest, we construct a model for the radial velocity distribution as a
function of $\phi_{1}$ (only DESI measurements are used in determining the
spline). We model the distribution as a Gaussian mixture where the center of
the Gaussian representing the GD-1 stream is allowed to vary with $\phi_{1}$.
The other Gaussian represents the background contamination. The variation of
the radial velocity of stream stars is represented by a spline with 6 knots.
The velocity dispersion of the stream is also allowed to vary and is
represented by a spline with 4 knots. The mean and the velocity dispersion of
the background contamination are fitted parameters. The specific likelihood
function for the radial velocity conditional on $\phi_{1}$ is given in
Equation 2 of Koposov et al. (2023). The model is implemented using the Stan
probabilistic programming language (Carpenter et al., 2017) and is sampled
using CMDstanpy package. The splines were implemented using the stan-splines
package using code very similar to that provided in the supplementary material
in Koposov et al. (2023). 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. Velocity prior for individual spline
knots is a normal distribution $\mathcal{N}(0,300)$; the prior on the
logarithm of the stream velocity dispersions at knots is
$\mathcal{N}(1.6,1.6)$ which is also uninformative. For the background
contamination the priors on the mean velocity, velocity dispersion and
velocity gradient are also normal distributions with large standard
deviations.

Figure 5: Top: DESI-observations in the GD-1 region. Middle measurement of
velocities at spline knots is shown by black points with error bars. The red
curves show the samples from the posterior distribution. The grey points show
the stars used to obtain the spline fit. Bottom: The velocity dispersion of
the GD-1 stream as measured by the spline model. The error bars show the
$1\sigma$ uncertainties in the velocity dispersion measurements derived from
the posterior samples at four points along the stream.

Figure 5 (top) shows the measured radial velocity $V_{GSR}$ vs. $\phi_{1}$ for
all DESI targets within $\pm 3^{\circ}$ from the stream track. The middle
panel shows the same points in grey and location of spline knots in black. The
uncertainties on the velocities are not visible as they are $\lesssim$ 15
km/s. Red lines show posterior samples from radial velocity curves.

The lower panel of Figure 5 shows the constraints on the velocity dispersion
along the stream in the region $-55^{\circ}\leq\phi_{1}\leq 0^{\circ}$. The
error-bars show the 1$\sigma$ uncertainties determined from the posterior
samples. The velocity dispersion measurements at different positions along the
stream obtained with Stan are given in Table 2. The velocity dispersion in the
$-40\leq\phi_{1}\leq-20$ part of the stream is around 3-4$\,\mathrm{km\
s}^{-1}$ measured with a small $0.6-1$ $\,\mathrm{km\ s}^{-1}$ uncertainty.
Higher velocity dispersion at $\phi=-46$ is consistent with the same velocity
dispersion within the errors. The measured dispersion at $\phi_{1}=-4$ is
somewhat lower (1 $\,\mathrm{km\ s}^{-1}$) but is still within 2 $\sigma$ from
the measurement at $\phi_{2}=-20$.

Table 2: Stream velocity dispersions
---
$\phi_{1}$ | $\sigma_{V_{GSR}}$ | $\epsilon_{\sigma_{V_{GSR}}}$
degrees | $\,\mathrm{km\ s}^{-1}$ | $\,\mathrm{km\ s}^{-1}$
-46.0 | 6.195 | 2.38
-32.0 | 4.24 | 0.99
-18.0 | 3.27 | 0.60
-4.0 | 0.97 | 0.88

While the Stan spline fit above is obtained by using only DESI velocity
measurements (since these have the smallest uncertainties), the spline in
Figure 4 is extended beyond the range of $\phi_{1}$ values modeled by the Stan
spline, by selecting knots by eye.

###  3.4 Metallicity distribution of stream stars

Figure 6: [Fe/H] distribution for stars in the GD-1 stream using the proper
motion cuts, parallax cut, CMD cut in the GD-1-region in previous figures
(grey). Black histogram shows stars with the additional constraint that they
have $\delta V_{GSR}<30$$\,\mathrm{km\ s}^{-1}$ from the stream track in
Figure 4(a). Only DESI measurements are used in this figure.

Figure 6 shows the distribution of DESI [Fe/H] measurements (from the RVS
pipeline) for all stars (grey) selected using the proper motion, parallax, and
CMD cuts using the selections described in Sections 3.1 and 3.2. Stars with
the additional constraint that they lie within 30 $\,\mathrm{km\ s}^{-1}$ from
the stream track in Figure 4 (Section 3.3) are shown in black. It is clear
from the black histogram in Figure 6 that only 2 stars have [Fe/H] $>-1.5$ and
one more has [Fe/H]$<-3.5$. We exclude these stars bringing the sample from
DESI to 143 (126 new, 17 improved). With these metallicity cuts 139 GD-1 stars
from SDSS+LAMOST remain, along with 43 stars from MMT (Bonaca et al., 2020).

The median metallicity of the entire DESI GD-1 sample as measured by the RVS
pipeline is -2.49$\pm 0.16$ (uncertainty determined by bootstrapping). Given
the uncertainty of 0.1-0.2 dex in the [Fe/H] values recovered by the RVS
pipeline, this is consistent with values from previous works: [Fe/H] =
-2.24$\pm 0.21$ (Malhan & Ibata, 2019), [Fe/H]=-2.3$\pm 0.1$ (Bonaca et al.,
2020), and [Fe/H] = -2.47 recently reported by Ibata et al. (2023).

The DESI EDR dataset contains spectra for 507 stars that are also in the
APOGEE DR17 catalog (Koposov et al., 2024). A detailed study of the random and
systematic errors on DESI stellar atmospheric parameters compared with the
corresponding APOGEE measurements is presented in section 5.1 of that paper.
Their analysis shows that the median DESI metallicity has a median offset
[Fe/H]${}_{\mathrm{DESI}}-$[Fe/H]${}_{\mathrm{APOGEE}}=0.06^{+0.13}_{-0.18}$
dex for the RVS pipeline and a median offset $0.03^{+0.18}_{-0.16}$ dex for
the SP pipeline. The results of the RVS and SP pipeline for very metal poor
stars have also been compared with follow up spectra with OSIRIS on Gran
Telescopio Canarias presented by Allende Prieto et al. (2023) who show that
the [Fe/H] values show excellent agreement with the OSIRIS values with a mean
difference of 0.04 dex and a standard deviation of 0.14 dex. Since the DESI
pipelines are still being improved and the systematic error likely depends on
the stellar population, we do not apply any corrections for systematic
uncertainties in [Fe/H] to the metallicities in this work.

Figure 7:  Stars in the GD-1 region with $|\delta\phi_{2}|<3^{\circ}$ from the
stream track with parallax, CMD, proper motion, metallicity
($-3.5<$[Fe/H]$<-1.5$)) colored by velocity relative to the velocity spline in
Figure 4(a). Points show stars observed by: DESI-EDR (circles), SDSS-
DR9+LAMOST-DR8 (diamonds), and MMT-Hectochelle (triangles). The orange box
shows the region we analyze for the cocoon (where most of the DESI
observations lie). Boundaries of DESI tiles (fields) are shown by magenta
curves. Figure 8: Distribution of DESI stars in $\delta V_{GSR}$ and
$\delta\phi_{2}$. The centrally concentrated component of the GD-1 stream and
the broader cocoon component are clearly visible.

##  4 Identification and characterization of the cocoon

As discussed in Section 1, an extended distribution called a “cocoon” has been
identified surrounding the GD-1 stream and a few other tidal streams. While
the origin of such cocoons is still uncertain (see Section 5), with the
increased number of radial velocities provided by DESI it is now possible to
identify and characterize the GD-1 cocoon in greater detail.

Our final selection (after imposing the parallax, proper-motion, CMD, RV and
metallicity cuts previously described) and with $|\delta\phi_{2}|\leq
3^{\circ}$ can be seen in Figure 7. This plot only shows observations with
full 6-D phase space and metallicity information. Malhan et al. (2019b)
identified the GD-1 cocoon primarily in the region
$-10^{\circ}<\phi_{1}<0^{\circ}$ based primarily on Gaia positional and proper
motion data. DESI spectroscopic measurements confidently confirm the existence
of the cocoon as a broad component in $\phi_{2}$ extending between
$-20^{\circ}<\phi_{1}<-5^{\circ}$.

We now analyze the phase space distribution of the stars in the region
identified by the orange box. In this work we restrict our analysis to this
range of $\phi_{1}$ since it is where the majority of the DESI spectroscopic
measurements exist. We do not claim that this is the full extent of the
cocoon. Spectroscopic observations at other positions along $\phi_{1}$ and
extending perpendicular to the stream in $\phi_{2}$ are required to
confidently assess its full extent of the cocoon and will be presented in
future work. The cocoon might well extend further down to
$\phi_{1}\sim-38^{\circ}$ but we have few DESI-EDR observations in this
region. This region also includes the ‘spur’ feature (seen as orange triangles
above the main stream track) and we choose to ignore it here since a specific
model has been postulated for the origin of this structure (Price-Whelan &
Bonaca, 2018; Bonaca et al., 2019).

Figure 8 shows a broader velocity range of stars in the orange box in Figure 7
to illustrate that the velocity range that we use ($\pm 30$$\,\mathrm{km\
s}^{-1}$) does indeed completely capture the cocoon region around the narrower
GD-1 stream. Since Figure 8 only includes data observed by DESI, the extent of
the cocoon in $\phi_{2}$ in this figure may be underestimated.

Figure 9 shows various phase space projections of stars in the orange box in
Figure 7 relative to the stream track (as identified by the 4 stream-track
splines previously defined in Table 1). In Figure 9 the left column shows (top
to bottom): a histogram of $\delta\phi_{2}$ and scatter plots of $\delta
V_{GSR}$, $\delta\mu_{\phi_{1}}$, $\delta\mu_{\phi_{2}}$ vs. $\delta\phi_{2}$.
In all plots in this column a centrally peaked distribution in
$\delta{\phi_{2}}$ is clearly distinguishable from a more extended
distribution in $\delta\phi_{2}$ corresponding to the cocoon. The plots in the
right hand column show (top to bottom): a histogram of $\delta V_{GSR}$ and
scatter plots of $\delta\mu_{\phi_{1}}$ and $\delta\mu_{\phi_{2}}$ vs. $\delta
V_{GSR}$. We also note that the plots show evidence for a greater spread in
$\delta V_{GSR}$ and $\delta\mu_{\phi_{1}}$ (two components of motion along
the stream) than in $\delta\mu_{\phi_{2}}$ (motion perpendicular to the stream
direction).

Figure 9: Stars from region in the orange box in Fig. 7 with the same
selections. DESI observations are in black, SDSS+LAMOST in blue (no MMT stars
lie in the cocoon region). In each scatter plot and histogram it is clear that
there are both narrow and broad components in $\delta\phi_{2}$, $\delta
V_{GSR}$ and $\delta\mu_{\phi_{1}}$. The overall width in
$\delta\mu_{\phi_{2}}$ (perpendicular to the stream track) is smaller than in
$\delta\mu_{\phi_{1}}$ (along the stream track).

###  4.1 Multi-Gaussian Analysis

We wish to separate the “thin stream” and “cocoon” components in order to
separately examine their phase space distributions. We use a Gaussian mixture
model to fit observed quantities with two univariate Gaussians (representing
the displacement of stars from the stream track in $\delta\phi_{2}$). One
Gaussian represents the “thin stream” component
$\mathcal{N}_{1}(\bar{x}_{1},\Sigma_{1})$ while the other represents the
cocoon $\mathcal{N}_{2}(\bar{x}_{2},\Sigma_{2})$, with an unknown mixing
fraction $Q$. We adapt the formalism of Equation 17 of Hogg et al. (2010) to
write the likelihood function as follows, where $x_{i}$ refer to individual
values of $\delta\phi_{2}$. We assume that positional uncertainty on all
$\delta_{\phi_{2}}$ values is zero since the Gaia positional uncertainty at
$G=21$ is 0.1 mas which is much smaller (by a factor $10^{7}$) than the spread
in $\delta\phi_{2}$. $Q$ is the fraction of stars in Gaussian component
$\mathcal{N}_{1}$ (“thin stream”). The likelihood function, assuming each
component is well described by a Gaussian, is given by,

$\displaystyle\mathcal{L}=$
$\displaystyle\prod_{i=1}^{N}Q\mathcal{N}_{1}(x_{i}|\bar{x}_{1},\Sigma_{1},\sigma_{xi})+(1-Q)\mathcal{N}_{2}(x_{i}|\bar{x}_{2},\Sigma_{2},\sigma_{xi}),$
(4) $\displaystyle\mathcal{L}=$
$\displaystyle\prod_{i=1}^{N}\frac{Q}{\sqrt{2\pi(\Sigma_{1}^{2}+\sigma_{xi}^{2})}}\exp{\left(-\frac{(x_{i}-\bar{x}_{1})^{2}}{2(\Sigma_{1}^{2}+\sigma_{xi}^{2})}\right)}$
$\displaystyle+\frac{(1-Q)}{\sqrt{2\pi(\Sigma_{2}^{2}+\sigma_{xi}^{2}})}\exp{\left(-\frac{(x_{i}-\bar{x}_{2})^{2}}{2(\Sigma_{2}^{2}+\sigma_{xi}^{2})}\right).}$

We estimate the five parameters in the above model ($\bar{x}_{1}$,
$\Sigma_{1}$, $\bar{x}_{2}$, $\Sigma_{2}$, $Q$) in two different ways. First,
we use direct maximization of the likelihood in Equation 4 with scipy.optimize
package. We determine the distribution of the best-fit parameters by
bootstrapping the observed sample 1000 times.

Second, we generate the posterior distribution of the model parameters using
the ensemble Markov Chain Monte Carlo sampling method emcee (Foreman-Mackey et
al., 2013). We adopt uniform flat priors for all the parameters
($\bar{x}_{1}$, $\Sigma_{1}^{2}$, $\bar{x}_{2}$, $\Sigma_{2}^{2}$, $Q$). The
posterior samples are generated with 64 walkers until the convergence of the
mean parameter values to 1%. Since the maximum autocorrelation length of the
chain $N_{\rm corr}\approx 80$ steps, we discard the first 1600 samples as
“burn-in” and thin the remaining chain by a factor $N_{\rm corr}$. It results
in about 6000 independent chain steps.

Both direct maximization of the likelihood and the MCMC method give nearly
identical result within the uncertainties. Both methods also yield similar
means and standard deviations of the two Gaussians in $\delta\phi_{2}$ and
hence we report the results of the Bayesian MCMC analysis. The resulting mean
($\overline{\delta\phi_{2}}$ and standard deviation
($\Sigma_{\delta\phi_{2}}$) of the Gaussians are reported in Table 3 (top 3
rows, of 2nd & 3rd columns). Uncertainties on all parameters in this table are
16th and 84th percentile credible intervals inferred from the posterior
probability distributions.

Table 3: Results of Univariate and Bivariate multi-Gaussian analysis
---
Parameter | Univariate | Bivariate
| Thin | Cocoon | Thin | Cocoon
$\overline{\delta\phi_{2}}$ [deg] | $-0.118\pm 0.017$ | $-0.149\pm 0.161$ | $-0.109\pm 0.017$ | $-0.165\pm 0.154$
$\Sigma_{\delta\phi_{2}}$ [deg] | $0.149\pm 0.013$ | $1.265\pm 0.122$ | $0.150\pm 0.014$ | $1.245\pm 0.117$
$\overline{\delta V_{GSR}}$ [$\,\mathrm{km\ s}^{-1}$] | … | … | $-1.37\pm 0.49$ | $1.23\pm 1.24$
$\Sigma_{\delta V_{GSR}}[\,\mathrm{km\ s}^{-1}]$ | … | … | $2.60\pm 0.65$ | $8.06\pm 0.98$
$Q_{\rm thin}$ | $0.63\pm 0.05$ | … | $0.60\pm 0.05$ | …

Figure 10: Top: Relative radial velocity $\delta V_{GSR}$ vs.
$\delta\phi_{2}$. Points show the thin and lower velocity dispersion stream
stars (black) and thick and higher velocity dispersion cocoon stars (red) as
determined by the bivariate double-Gaussian model. The vertical lines are at
$\overline{\delta\phi}_{2,thin}\pm\Sigma_{\delta\phi_{2,thin}}$ and include
the majority of the thin stream stars identified as members of the thin
component by both position and velocity. Bottom: Kernel density estimates of
in $\delta\phi_{2}$ for all stars in the orange box in Figure 6 with Gaussians
representing the thin stream (black curve) and cocoon (red curve) as estimated
by the bivariate analysis. Figure 11: Relative radial velocity and both
components of proper motions of the stars from stream (left column) and cocoon
(right column) using membership probabilities from the bivariate Gaussian
model to define thin and cocoon components.

We also fit the $\phi_{2}$ distribution and the $V_{GSR}$ distribution jointly
with two bivariate-Gaussian distributions described by the likelihood function

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\prod_{i=1}^{N}Q\mathcal{N}_{1}(x_{i},y_{i},|\bar{x}_{1},\bar{y}_{1},\Sigma_{x1},\Sigma_{y1},\sigma_{xi},\sigma_{yi})+$
(5)
$\displaystyle(1-Q)\mathcal{N}_{2}(x_{i},y_{i},|\bar{x}_{2},\bar{y}_{2},\Sigma_{x2},\Sigma_{y2},\sigma_{xi},\sigma_{yi}),$

where both $\mathcal{N}_{1}$ (thin stream) and $\mathcal{N}_{2}$ (cocoon) are
represented by 2D (bivariate) Gaussians $\mathcal{N}(\delta\phi_{2},\delta
V_{GSR})$. The results of the Bayesian MCMC analysis of the bivariate Gaussian
model (9 parameters and their uncertainties) are also given in Table 3.

Both the univariate and bivariate analyses give nearly identical widths in
$\delta\phi_{2}$ for the ‘thin stream’ and cocoon. The standard deviation of
the Gaussian fit to the thin component $\Sigma_{\delta\phi_{2},thin}=(0.15\pm
0.014)^{\circ}$ implies a FWHM of the “thin stream” component in $\phi_{2}$ to
be 2.355$\times\sigma_{\phi_{2},thin}=0.353^{\circ}\approx 55$ pc at the
distance of 9 kpc) while the cocoon component is much broader with standard
deviation $\Sigma_{\delta\phi_{2},\mathrm{cocoon}}=(1.245\pm 0.117)^{\circ}$
or FWHM $\approx 2.931^{\circ}\approx 460$ pc. For the distribution in
$V_{GSR}$ the bivariate analysis gives velocity dispersions
$\Sigma_{V_{GSR},\mathrm{thin}}=(2.60\pm 0.65)$ $\,\mathrm{km\ s}^{-1}$ and
$\Sigma_{V_{GSR},\mathrm{cocoon}}=(8.06\pm 0.98)$ $\,\mathrm{km\ s}^{-1}$ for
the thin stream and cocoon components respectively, with the fraction of stars
in the thin stream $Q=0.60\pm 0.05$.

The standard deviations in $\delta\phi_{2}$ of the thin stream and cocoon
obtained both by only modeling $\delta\phi_{2}$ versus by jointly modeling
$\delta\phi_{2}$ and $\delta V_{GSR}$ yields nearly identical results for the
angular widths of the thin stream and cocoon. Both methods also give $Q\sim
0.6$ for the fraction of stars in the ‘thin stream’ component.

Figure 10 (upper panel) shows the stars belonging to the “thin stream” with
probability $p_{thin}\geq 0.5$ as determined by the bivariate analysis
considering both position and velocity $\delta V_{GSR}$. The vertical lines
are located at $\overline{\delta\phi}_{2,thin}\pm\Sigma_{\delta\phi_{2,thin}}$
from the bivariate analysis in Table 3 and include the majority of points
determined to be members of the “thin” component by both univariate and
bivariate analyses. The vertical lines are not used in subsequent analysis and
only meant to guide the eye. Figure 10 (lower panel) shows the kernel density
estimate of all the stars in the orange box in Figure 7. Overplotted as black
and red curves are the Gaussians in $\delta\phi_{2}$ for the thin stream and
cocoon respectively (as determined by the bivariate Gaussian analysis).

Figure 11 shows various phase-space projections of the thin stream (left) and
the broader cocoon component (right) with membership defined by the individual
probabilities obtained from the bivariate analysis. It is clear that the
cocoon has a much greater spread in $\delta\mu_{\phi_{1}}$ and $\delta
V_{GSR}$ than in $\delta\mu_{\phi_{2}}$. Although the “thin stream” component
defined in this way has a very narrow central component in
$\delta\mu_{\phi_{1}}$, some stars show significant velocity and proper motion
deviations. The “thin stream” component also appears to have both narrow and
broad components in $\delta V_{GSR}$. Both the thin stream and cocoon
components have narrow distributions in $\delta\mu_{\phi_{2}}$, but the cocoon
has broad distributions in both $\delta V_{GSR},\delta\mu_{\phi_{1}}$.

###  4.2 Accounting for contamination by background stars

In the above analysis we have not explicitly modeled the background
population. Adding a third bivariate Gaussian component to represent
contamination by a background population adds 5 additional parameters. Given
the small number of stars in total in this region of the stream we did not
attempt such a model. Instead we use the fact that both of the previous
analyses yielded consistent values for the width of the thin stream. We use
the probabilities of membership to the ‘thin stream’ and cocoon based on the
univariate analysis in Section 4.1 to define each component. We then model
each component’s velocity distribution with a univariate Gaussian component
for the velocity of the thin stream (cocoon). We represent the background by a
much broader Gaussian in velocity truncated at $\delta V_{GSR}=\pm
30$$\,\mathrm{km\ s}^{-1}$ (the velocity cut shown by the dashed lines in Fig.
4a). The mean and velocity dispersion of the background contamination are
fitted parameters. The likelihood function is similar to that given in
Equation 4.

As one would expect if some fraction of the higher velocity stars in both the
‘thin stream’ and cocoon are classified as belonging to the background, the
velocity dispersions of these two components would decrease.

This model yields a velocity dispersion of =$2.24_{-0.41}^{+0.48}$
$\,\mathrm{km\ s}^{-1}$ with the fraction of stars assigned to the background
found to be $0.15\pm 0.04$ This velocity dispersion for the thin stream
component is consistent (within the uncertainies) to the dispersion obtained
from the bivariate multi-Gaussian analysis in the previous section.

When applied to the members of the cocoon component we find that the cocoon
has a smaller velocity dispersion of $4.94^{+0.78}_{-0.66}$ $\,\mathrm{km\
s}^{-1}$ and the background fraction in the cocoon is estimated to be 12%.

To assess whether this background fraction is reasonable we counted the number
of stars (after applying the same cuts used to select stars in the orange box
in Fig. 7) in three separate off-stream regions with the same size as the
orange box. The number of stars in these 3 regions compared to the number of
cocoon stars varied from 4-8% suggesting that the number of background stars
in the cocoon could be somewhat lower than in the above estimate of 12%.

The analyses in this subsection and the previous one both imply that the
physically thin stream has a velocity dispersion of $\sim
2.2-2.6$$\,\mathrm{km\ s}^{-1}$, while the cocoon has a velocity dispersion
between 5-8$\,\mathrm{km\ s}^{-1}$. Additional data from future DESI data
releases will enable us to determine if the cocoon is an extended structure
found along the entire stream or if it is more localized. It will also enable
us to derive a more robust estimate of the velocity dispersion of the cocoon.

Our measurement of the dispersion of the ‘thin’ stream component
($\sigma_{V_{GSR},\mathrm{thin}}\sim 2$ $\,\mathrm{km\ s}^{-1}$) is consistent
with previous estimates of both other GC streams and the GD-1 stream itself
(Gialluca et al., 2021). In contrast, the larger velocity dispersion of the
broader cocoon component of $\sigma_{V_{GSR},\mathrm{cocoon}}\sim 5-8$
$\,\mathrm{km\ s}^{-1}$ is comparable to the velocity dispersion of classical
dwarf spheroidal galaxies like Carina, Draco, Fornax, Leo I etc. (Simon,
2019). This measurement is also comparable to the velocity dispersion of other
dwarf galaxy streams (e.g., Li et al., 2022) which have velocity dispersions
that range from 4-20$\,\mathrm{km\ s}^{-1}$ with a mean value of around
11$\,\mathrm{km\ s}^{-1}$ and are significantly larger than the dispersions of
known GC streams (e.g., Li et al., 2022).

###  4.3 Metallicity of thin stream and cocoon components

Another way to assess if the cocoon velocity dispersion in Section 4.1 is
inflated by background/foreground halo stars is to examine the metallicity
distributions of stars in the two kinematically and spatially identified
components.

Figure 12: Kernel density estimates (with an Epanechnikov kernel) of [Fe/H] as
reported by the DESI-RVS pipeline for the ‘thin stream’ component (black/grey)
and the cocoon component (red). The average random error on the [Fe/H]
measurements is $\sim 0.2$ dex and is also the band-width of the kernel
density estimate. Overplotted are the best-fit Gaussians (see text for
details).

Figure 12 shows the [Fe/H] distribution functions for the thin stream
component (grey/black) and the cocoon component (red), where the membership
probability from the bivariate analysis was used to assign memberships. It is
clear that the distribution of [Fe/H] for cocoon stars overlaps significantly
with the distribution for thin stream stars.

Once again assuming that the metallicity distribution for each component is
described by a Gaussian we use a Bayesian MCMC analysis (32 walkers, 1000
steps) to determine their intrinsic means and intrinsic standard deviations
while accounting for the uncertainties in the measurements. The mean
metallicity of the ‘thin stream’ component is found to be
$\left<\mathrm{[Fe/H]_{thin}}\right>=-2.54\pm 0.04$dex and intrinsic
metallicity dispersion is $\sigma_{\mathrm{[Fe/H],thin}}=0.27\pm 0.03$dex.

The cocoon component has a mean
$\left<\mathrm{[Fe/H]_{\mathrm{cocoon}}}\right>=-2.45\pm 0.06$dex with
dispersion $\sigma_{\mathrm{[Fe/H],\mathrm{cocoon}}}=0.30^{+0.05}_{-0.04}$dex.
The means of both these distributions are consistent with the value of [Fe/H]
= $-$2.47 dex recently reported by Ibata et al. (2023) for the entire GD1
stream. They are also significantly lower than the stellar halo metallicity of
$\sim-1.5$ dex. A Kolomogorov-Smirnov 2-sided test yields $p=0.15$, implying
that we are unable to reject the null hypothesis that the two distributions
are drawn from the same parent distribution.

A comparison with other Milky Way tidal streams (see Table 1 of, Li et al.,
2022) shows that the inferred instrinsic metallicity dispersion of GD-1 is
significantly higher than typical GC streams which have metallicity
dispersions of 0.03-0.05dex. This rather high inferred intrinsic metallicity
dispersion of GD-1 can be completely explained by the current systematic floor
in the uncertainty of [Fe/H] measurements resulting from the RV pipeline.
Cooper et al. (2023) showed that [Fe/H] measurements with DESI for several GCs
(M13, M92, NGC 5053) showed significant spreads of 0.3 dex or greater and
often had mean metallicities that were systematically lower by 0.2 to 0.4 dex
relative to literature values.

We caution that the DESI EDR metallicity values should be treated as
preliminary until ongoing improvements to the spectroscopic pipelines are
completed. Some of the improvements in progress including machine learning
based methods are expected to yield more accurate metallicities and abundances
from the low resolution, low S/N DESI spectra (J. Han, & DESI Collaboration,
in prep, M. Haiger, & DESI Collaboration, in prep).

##  5 Discussion and conclusions

###  5.1 Summary of results

We have spectroscopically identified 126 new members of the GD-1 tidal stream
and improved RV measurements for 17 more using observations with Dark Energy
Spectroscopic Instrument (DESI Collaboration et al., 2016, 2022).

We identify stars associated with the GD-1 stream by selecting in position
$|\delta\phi_{2}|<3^{\circ}$, and proper motion $|\delta\mu_{\phi_{1}}|<2$mas
yr-1, $|\delta\mu_{\phi_{2}}|<2$mas yr-1 relative to the stream track. Next we
apply a selection using a PARSEC isochrone with [M/H] = $-$2 and age = 12 Gyr
and selecting stars with $|\delta(g-r)|<0.1$ from this isochrone. After this
selection we identify and fit a spline to the stream-track in
$\phi_{1}-V_{GSR}$ coordinates using a Gaussian mixture modelling approach,
accounting for contamination by a background population, following the
approach described in Koposov et al. (2023). Rejecting two stars with
[Fe/H]$>-1.5$ and one more with [Fe/H]$<-3.5$ brings the total number of DESI
observations of GD-1 members to 143 of which 126 are new members and for the
remaining DESI has decreased the radial velocity uncertainty. Slightly
relaxing the cuts we made in proper motion, color-magnitude, and the choice of
the stellar isochrone used, can all increase the numbers of stars by 10-30
stars depending on the specific choices made. We consider our selection to be
conservative.

With these new data we spectroscopically identify a broad cocoon structure
surrounding a thinner and kinematically colder stream in the region
$-25^{\circ}\leq\phi_{1}\leq-5^{\circ}$ of the GD-1 tidal stream (where the
majority of our observations were obtained). This cocoon was first identified
by Malhan et al. (2019b) using primarily Gaia positions and proper motions and
was recently confirmed by Ibata et al. (2023).

Using a bivariate Bayesian mixture modeling algorithm and assuming that the
thin stream and cocoon are described by Gaussians, we determine their angular
widths to be $\Sigma_{\delta\phi_{2},\mathrm{thin}}=0.150\pm 0.014^{\circ}$
and $\Sigma_{\delta\phi_{2},\mathrm{cocoon}}=1.245\pm 0.117^{\circ}$. Since
our DESI fields are restricted to $|\delta_{\phi_{2}}|\sim\pm 3^{\circ}$ from
the stream track this could be an under-estimate of the width of the cocoon.
The velocity dispersions of these two components are
$\Sigma_{V_{GSR},\mathrm{thin}}=2.60\pm 0.65$$\,\mathrm{km\ s}^{-1}$ and
$\Sigma_{V_{GSR},\mathrm{cocoon}}=8.06\pm 0.98$$\,\mathrm{km\ s}^{-1}$. An
analysis of the velocity dispersions of the thin stream and cocoon components
defined by membership probability using $\delta\phi_{2}$ only, but allowing
for contamination by a background population gives a nearly identical velocity
dispersion for the thin stream component but a somewhat lower velocity
dispersion of $\sim 5$$\,\mathrm{km\ s}^{-1}$ for the cocoon component. A
bivariate analysis also including a background would benefit from more data
and is deferred to a future work.

DESI spectroscopy yields mean [Fe/H] values of
$\left<\mathrm{[Fe/H]_{thin}}\right>$ $=-2.50\pm 0.03$dex and
$\left<\mathrm{[Fe/H]_{\mathrm{cocoon}}}\right>$ $=-2.48\pm 0.03$dex
respectively, with intrinsic metallicity dispersions
$\sigma_{\mathrm{[Fe/H],thin}}=0.27\pm 0.02$dex and
$\sigma_{\mathrm{[Fe/H],\mathrm{cocoon}}}=0.29\pm 0.03$dex respectively. Our
measurements of the mean metallicity are consistent with previous estimates of
[Fe/H] for the whole GD-1 stream and also imply that the cocoon and the thin
stream probably originated from the same progenitor. The measured metallicity
dispersions are likely overestimated and will be improved.

###  5.2 Origin of the “cocoon”

While its low metallicity, retrograde orbit and the thinness and coldness of
the most visible part of the stream clearly imply that the progenitor of the
GD-1 stream was an accreted GC, the origin of the cocoon is currently
uncertain. At least four possibilities exist.

  1. 1.

The cocoon could comprise GC stars that were tidally stripped from the GC
while it was still inside its parent dwarf galaxy (pre-accretion tidal
stripping);

  2. 2.

it could be a tidal stream from a very metal poor dwarf galaxy, while the thin
component could be the stream from an associated nuclear star cluster or GC;

  3. 3.

a thin tidal stream generated from the GC progenitor could have been heated by
multiple interactions with dark subhalos in the MW;

  4. 4.

a thin tidal stream generated from the GC progenitor could have been perturbed
and heated by a massive Sagittarius dwarf galaxy.

We now describe each of these possible origins in turn and discuss how one
might distinguish between them with future spectroscopic data.

#### Cocoon stars originated from pre-accretion tidal stripping of the
progenitor GC:

Simulations show that star clusters can experience the strongest tidal force
in the first few hundred Myr after formation (Meng & Gnedin, 2022). Once the
parent dwarf galaxy is accreted onto the Milky Way and is tidally disrupted,
the pre-accretion tidal debris forms the cocoon while the intact GC is also
liberated from the parent dwarf galaxy and continues to be tidally stripped
forming the thinner GC stream (Carlberg, 2018; Malhan et al., 2019b, 2020;
Carlberg & Agler, 2023).

The structure and kinematics of cocoons so formed may be important probes of
the central dark matter distribution of the parent dwarf galaxy. (Malhan et
al., 2020) carried out a suite of over 100 controlled N-body simulations of
dwarf galaxies and their GCs falling into a Milky Way-like potential. They
showed that accreted GCs born in dwarf galaxies with cored dark matter halos
exhibit fewer morphological peculiarities (e.g., gaps and off-track stream
structures like the spur and blob) than their counterparts born in dwarf
galaxies with cuspy dark matter halos. They argue that for a parent halo of a
given mass, the physical width, radial and tangential velocity dispersions and
spread in angular momentum averaged along an accreted GC stream can serve as
diagnostics of whether the parent host galaxy had a cusp or a core. Malhan et
al. (2020) showed that if considering dwarf galaxies of halo mass
$10^{9}$$\,M_{\odot}$, the gentler tidal fields in a cored dark matter halo
would result in tidal stream with a velocity dispersion of $\sim
2-5$$\,\mathrm{km\ s}^{-1}$, while a dwarf galaxy with a cuspy halo would
produce a stream velocity dispersion of order 8-12$\,\mathrm{km\ s}^{-1}$.

In this pre-accretion tidal stripping scenario a cocoon velocity dispersion of
$\sim 8$$\,\mathrm{km\ s}^{-1}$ would imply that GD-1’s parent dwarf galaxy
had a cuspy dark matter halo while $\sim 5$$\,\mathrm{km\ s}^{-1}$ would imply
the parent dwarf had a cored halo. Additional observations along a larger
stretch of the stream are required to draw definitive conclusions since this
formation scenario predicts that cocoon are not restricted to short segments
of the stream.

#### Cocoon stars originated from the parent dwarf galaxy:

A second possibility is that the cocoon comprises the stars of GD-1’s parent
dwarf galaxy. The velocity dispersion of the cocoon as well as the metallicity
are reasonably consistent with those of other Local Group dwarf galaxies. If
the cocoon is the debris of the parent dwarf galaxy, the metallicity of the
cocoon would imply that the parent had a stellar mass of $\sim
10^{4}-10^{6}$$\,M_{\odot}$(Kirby et al., 2013). Using the stellar mass - halo
mass relation for low mass halos (Zaritsky & Behroozi, 2023) this would imply
a parent dwarf galaxy of halo mass $\sim 2\times 10^{9}$$\,M_{\odot}$.

Dwarf galaxies tend to show greater metallicity dispersions than GCs but the
metallicity distributions of the thin stream and cocoon components are
currently statistically indistinguishable, although improved [Fe/H]
measurements are needed to make a definitive statement. Nuclear star clusters
do show large metallicity dispersions (e.g., Alfaro-Cuello et al., 2019;
Nogueras-Lara, 2022), but all known nuclear star clusters are more metal rich
than GD-1.

GD-1’s low metallicity and retrograde orbit are strong indicators that
progenitor was accreted from a dwarf galaxy. Previous work based on the
location of GD-1 in action space relative to known GCs and halo stars suggests
a possible association with debris from accreted dwarf galaxy(ies) such as
‘Sequoia/Arjuna/ I’itoi’ (Myeong et al., 2019; Bonaca et al., 2021; Malhan et
al., 2022b). Accreted GCs which are about 60% of the GC population (Choksi et
al., 2018; Choksi & Gnedin, 2019; Creasey et al., 2019) tend to be on more
radial orbits and tend to be found at larger radii (Chen & Gnedin, 2022, 2023;
Belokurov & Kravtsov, 2024). In contrast, in situ GCs are more metal-rich and
lie deeper in the Galactic potential (at more negative energy). GD-1 is
peculiar because it is quite deep in the potential but rather than being on a
radial orbit it is on a retrograde and nearly tangential (high angular
momentum) orbit. A possible explanation is that it was brought in by a fairly
heavy dwarf galaxy that experienced significant dynamical friction (Binney &
Tremaine, 2008).

Follow-up high-resolution spectroscopy to obtain abundance data will be
particularly useful in determining whether the cocoon has abundance signatures
characteristic of GCs or dwarf galaxies. For example a sodium-oxygen anti-
correlation is seen in high-resolution spectra of most GCs (e.g., Marino et
al., 2011) but not in dwarf galaxies.

#### Cocoon produced by dark matter subhalo impacts with a thin GC stream:

The third possible origin for the cocoon is heating by impacts from multiple
dark subhalos with masses in the range $10^{7}-10^{8}$$\,M_{\odot}$ (e.g.,
Johnston et al., 2002; Ibata et al., 2002; Siegal-Gaskins & Valluri, 2008;
Carlberg, 2012, 2013; Erkal & Belokurov, 2015; Sanders et al., 2016; Banik et
al., 2018; Bonaca et al., 2019; Banik et al., 2021; Carlberg & Agler, 2023).
Simulations provide insight into the origins of the broad cocoon-like
components produced by subhalo impacts. Encounters with small impact
parameters between thin GC streams and dark subhalos cause the decrease of the
angular momenta of stream stars ahead of the point of closest approach, since
stars are pulled towards the perturber, and vice-versa. Stars with lower
angular momenta than the progenitor spread to smaller galactocentric radii
while stars with higher angular momenta spread to larger radii. The
accumulation of sub-halo encounters with time is expected to lead to an
increase in velocity dispersion (Ngan et al., 2015, 2016) along the stream and
may also create gaps which orbital shearing will disperse over time. In
addition, orbital precession in the non-spherical potential of a Milky Way-
like system causes orbits to diverge on a time-scale of several Gyr (Ibata et
al., 2002; Siegal-Gaskins & Valluri, 2008). Since the effects of sub-halo
encounters are initially fairly local, only old GC streams (from star clusters
accreted long ago) will have had time to experience subhalo heating over
significant stretches of the stream. Recent models suggest that the GD-1
stream is between 3-5.5 Gyr old (Gialluca et al., 2021), suggesting it might
not be old enough to have been heated by subhalos along its entire extent. On-
going and future observations of the GD-1 stream with DESI will enable us to
determine whether the cocoon is a local or a much more extended feature.

#### Perturbations from a massive Sagittarius dwarf:

de Boer et al. (2020) proposed that the spur in GD-1 could have been caused by
interaction of the Sagittarius dwarf spheroidal galaxy (Sgr) with the GD-1
stream. Dillamore et al. (2022a) showed that if the Sgr dwarf was originally
more massive than $4\times 10^{10}$$\,M_{\odot}$, gravitational torques from
it would affect tidal streams (like Pal-5 and GD-1) that lie within the
orbital pericenter radius of Sgr. Such torques can cause asymmetries between
the leading and trailing arms due to the folding back of one of the arms as
well as the appearance of off-stream structures like the GD-1 “spur” and
“blob” (Price-Whelan & Bonaca, 2018) (the “blob” overlaps with the cocoon we
characterized in previous subsections). Perturbations from a massive Sgr may
even produce structures that appear like a separate stream on the sky but with
kinematics similar to the main stream. Dillamore et al. (2022a) show that if
Sgr has affected the GD-1 stream it could have also produced a thin off-set
stream similar to Kshir, a stream intersecting GD-1 discovered by (Malhan et
al., 2019a) .

In practice more than one of the above processes could have affected the GD-1
stream to varying degrees: the cocoon could have contributions from both pre-
accretion tidal heating, stars from the parent dwarf galaxy and heating from
sub-halo encounters. Additional observations, as well as modeling, will help
to untangle the contributions from various processes.

Recently Ibata et al. (2023) also identified a two-component structure for the
GD-1 stream. They used the VLT/UVES spectrograph to follow up potential GD-1
candidates detected in the Gaia EDR3 and DR3 catalogs. The VLT/UVES data have
velocity precision of 1$\,\mathrm{km\ s}^{-1}$. They also obtained
observations with the IDS long-slit spectrograph on the 2.5m Isaac Newton
Telescope which have a lower velocity precision ($\sim 10$$\,\mathrm{km\
s}^{-1}$). Using 323 stars with line of sight velocity measurements they find
a two component fit with velocity dispersions for the “thin stream” to be
$\sigma=7.4\pm 1.1$$\,\mathrm{km\ s}^{-1}$ while the cocoon has a velocity
dispersion $\sigma=29.1\pm 6.1$$\,\mathrm{km\ s}^{-1}$ with a fraction of
$0.73\pm 0.09$ in the thin component.

Our analysis in section 4 is restricted to a narrow range of $\phi_{1}$ values
and implies a much smaller velocity dispersion of 5-8$\,\mathrm{km\ s}^{-1}$
than found by Ibata et al. (2023). Their analysis is applied to the entire
stream ($-100^{\circ}<\phi_{1}<20^{\circ}$) while our analysis above is
restricted to $-30^{\circ}<\phi_{1}<0^{\circ}$. Although we have few DESI
observations at $\phi_{1}<-50^{\circ}$, when we use data from DESI,
SDSS+LAMOST, MMT (all the data within the dashed red curves in Figure 4) to
estimate the width of the thin stream and cocoon we find that the velocity
dispersions of the thin component and cocoon are $1.89\pm 0.38$$\,\mathrm{km\
s}^{-1}$ and $7.68\pm 0.72$$\,\mathrm{km\ s}^{-1}$ respectively with a thin
stream fraction $Q=0.58$. Although these values are probably dominated by our
measurements in the cocoon area, they are significantly smaller than the
estimates from Ibata et al. (2023). While our estimated velocity dispersion
for the thin stream component is more consistent with previous work on the
GD-1 stream’s thin component as well as other GC streams. Our estimated
velocity dispersion for the cocoon is also more consistent with simulation
predictions (Malhan et al., 2020; Carlberg & Agler, 2023).

At the present time it is not possible to draw definitive conclusions as to
the origin of the cocoon component of the GD-1 stream. The more important
conclusion is that the GD-1 stream does have narrow and wide components
approximately as expected in Milky Way simulations with an accreting sub-halo
population. These initial data suggest that definitive conclusions will be
possible as increased numbers of stars in individual streams become available
and more streams are studied.

###  5.3 Other factors affecting the structure of the GD-1 stream

In addition to the mechanisms discussed in the previous section, there have
been recent studies of at least a few other factors that can affect the
structure and kinematics of tidal streams like GD-1.

It is thought that the massive merger in the Milky Way’s history referred to
as the Gaia-Enceladus-Sausage event (Belokurov et al., 2018; Helmi et al.,
2018) could have resulted in the disk of the Milky Way becoming misaligned
with the resulting triaxial dark matter halo (Dillamore et al., 2022b; Han et
al., 2022, 2023) or otherwise cause torques on the stellar disk which cause it
to tilt as it realigns with the halo. Simulations show that the tilting rate
of disks following mergers can have angular speeds of $\sim
10-20$$\,\mathrm{km\ s}^{-1}$kpc-1 (e.g., Earp et al., 2017; Dodge et al.,
2023). Recently Nibauer et al. (2023) showed that the gravitational torque
from such a tilting disk can cause both narrowing and broadening (depending on
the magnitude and direction of tilting) of tidal streams with pericenters
within $\sim 20\,$kpc. Tilting can also cause changes in the line-of-sight
width $\delta\phi_{2}$ as well as the width of the radial velocity
distribution. They show that several of the features of the Pal-5 stream are
consistent with a disk tilting rate of 15$\,\mathrm{km\ s}^{-1}$kpc-1. While
these authors do not discuss the GD-1 stream it could also be affected by disk
tilting since its pericenter radius is $<20$ kpc.

Despite the large number of possible origins for the cocoon described in this
and the previous subsection, detailed models and additional spectroscopic data
may enable us to determine its origin.

###  5.4 GD-1 data expected from future DESI data releases

Future DESI data releases will include many more GD-1 members. The first year
of the DESI main program (Year 1) has already obtained spectra for $\sim$ 350
GD-1 members from a total of $\sim 1300$ targets observed in the stream area
(2000 square degrees). Targets were selected by CMD and proper motion cuts and
most were observed as part of secondary and tertiary observing programs
designed for future survey planning purposes. The regions with confirmed
members span $-75^{\circ}<\phi_{1}<15^{\circ}$. The sample of confirmed
members include stars down to $r=21$ (the velocity uncertainty at the faint
end is $\sim 15$$\,\mathrm{km\ s}^{-1}$). Current plans for the future include
observing up to ninety 3∘ diameter fields in the GD-1 stream region. It is
anticipated that by the end of the 5-year program DESI will obtain spectra for
$\sim 3000$ GD-1 stream members, making it one of the most thoroughly surveyed
GC streams.

MV gratefully acknowledges financial support from NASA-ATP award
80NSSC20K0509. SK acknowledges support from the Science & Technology
Facilities Council (STFC) grant ST/Y001001/1. T.S.L. acknowledges financial
support from Natural Sciences and Engineering Research Council of Canada
(NSERC) through grant RGPIN-2022-04794. APC acknowledges support from a Taiwan
Ministry of Education Yushan Fellowship and the Taiwan National Science and
Technology Council (NSTC) grants 109-2112-M-007-011-MY3 and
112-2112-M-007-017-MY3. This material is based upon work supported by the U.S.
Department of Energy (DOE), Office of Science, Office of High-Energy Physics,
under Contract No. DE–AC02–05CH11231, and by the National Energy Research
Scientific Computing Center, a DOE Office of Science User Facility under the
same contract. Additional support for DESI was provided by the U.S. National
Science Foundation (NSF), Division of Astronomical Sciences under Contract No.
AST-0950945 to the NSF’s National Optical-Infrared Astronomy Research
Laboratory; the Science and Technology Facilities Council of the United
Kingdom; the Gordon and Betty Moore Foundation; the Heising-Simons Foundation;
the French Alternative Energies and Atomic Energy Commission (CEA); the
National Council of Science and Technology of Mexico (CONACYT); the Ministry
of Science and Innovation of Spain (MICINN), and by the DESI Member
Institutions: https://www.desi.lbl.gov/collaborating-institutions. Any
opinions, findings, and conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect the views
of the U. S. National Science Foundation, the U. S. Department of Energy, or
any of the listed funding agencies. The DESI Legacy Imaging Surveys consist of
three individual and complementary projects: the Dark Energy Camera Legacy
Survey (DECaLS), the Beijing-Arizona Sky Survey (BASS), and the Mayall z-band
Legacy Survey (MzLS). DECaLS, BASS and MzLS together include data obtained,
respectively, at the Blanco telescope, Cerro Tololo Inter-American
Observatory, NSF NOIRLab; the Bok telescope, Steward Observatory, University
of Arizona; and the Mayall telescope, Kitt Peak National Observatory, NOIRLab.
NOIRLab is operated by the Association of Universities for Research in
Astronomy (AURA) under a cooperative agreement with the National Science
Foundation. Pipeline processing and analyses of the data were supported by
NOIRLab and the Lawrence Berkeley National Laboratory (LBNL). The Legacy
Surveys also use data products from the Near-Earth Object Wide-field Infrared
Survey Explorer (NEOWISE), a project of the Jet Propulsion
Laboratory/California Institute of Technology, funded by the National
Aeronautics and Space Administration. Legacy Surveys was supported by: the
Director, Office of Science, Office of High Energy Physics of the U.S.
Department of Energy; the National Energy Research Scientific Computing
Center, a DOE Office of Science User Facility; the U.S. National Science
Foundation, Division of Astronomical Sciences; the National Astronomical
Observatories of China, the Chinese Academy of Sciences and the Chinese
National Natural Science Foundation. LBNL is managed by the Regents of the
University of California under contract to the U.S. Department of Energy. This
work has made use of data from the European Space Agency (ESA) mission
Gaia777https://www.cosmos.esa.int/gaia, processed by the Gaia Data Processing
and Analysis Consortium
(DPAC888https://www.cosmos.esa.int/web/gaia/dpac/consortium).This work has
made use of data from the Sloan Digital Sky Survey. Funding for SDSS-III has
been provided by the Alfred P. Sloan Foundation, the Participating
Institutions, the National Science Foundation, and the U.S. Department of
Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. This
work also utilized data from the LAMOST survey. This paper made use of the
Whole Sky Database (wsdb) created and maintained by Sergey Koposov at the
Institute of Astronomy, Cambridge, with financial support from STFC and the
European Research Council (ERC). The authors are honored to be permitted to
conduct scientific research on Iolkam Du’ag (Kitt Peak), a mountain with
particular significance to the Tohono O’odham Nation. 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: Data behind all figures in this paper are available as FITS
files at: https://zenodo.org/records/11638330 DOI: 10.5281/zenodo.12535192.

##  Appendix A Measuring the distance stream track

To calibrate the distance gradient along the stream, we use the proper motion
and parallax selection described in Section 3 and then focus on the subgiant
branch stars as those form a narrow sequence on the CMD (see Figure 3).
Specifically we select stars with $0.25<g-r<0.4$ and assume that the subgiant
branch has an approximately constant $r+4(g-r)$. The figure 13 shows the
$r+4(g-r)-4.5$ vs. $\phi_{1}$ for likely GD-1 subgiant stars. The red curve
shows our fiducial distance modulus vs. $\phi_{1}$ relation.

Figure 13: Left: The stream distance gradient from the subgiant branch stars.
The greyscale shows the distribution of the $r+4(g-r)-4.5$ (a proxy for
distance modulus) as a function of $\phi_{1}$ for stars with $.25<g-r<.4$. The
subgiant branch stars form a narrow sequence in the figure. The red curve
shows the adopted distance gradient model (shifted down by 0.5 magnitudes for
clarity). Right: Various relations for the distance to the GD-1 stream as a
function of $\phi_{1}$ from Price-Whelan & Bonaca (2018) (blue dot-dashed), de
Boer et al. (2018) (green dashed), Li et al. (2018) (black solid) and this
work (red solid).

##  Appendix B Table of DESI-EDR measurements in the GD-1 area

Table 4 gives DESI-EDR measurements of $V_{GSR}$ and [Fe/H] and their
1$\sigma$ measurement uncertainities, $V_{err}$ and [Fe/H]err are from the RVS
pipeline. The $V_{GSR}$ measurements have been corrected by 0.93$\,\mathrm{km\
s}^{-1}$, the systematic error on RV measurements reported in (Koposov et al.,
2024). No corrections have been applied to [Fe/H] measurements. TARGET_ID is
the DESI target ID, RA and Dec values are from Gaia-DR3.

Table 4: DESI-EDR Observations of GD-1 region (data in Fig 4a)
---
TARGET-ID | RA | Dec | $\phi_{1}$ | $\phi_{2}$ | $V_{GSR}$ | $V_{err}$ | [Fe/H] | [Fe/H]err
| [deg] | [deg] | [deg] | [deg] | [$\,\mathrm{km\ s}^{-1}$] | [$\,\mathrm{km\ s}^{-1}$] | [dex] | [dex]
39633374… | 192.76… | 58.91… | -3.70… | 0.49… | 125.44… | 1.91… | -1.92… | 0.05
39633362… | 191.93… | 57.92… | -4.44… | -0.29… | 113.64… | 1.49… | -1.50… | 0.036
39633377… | 191.69… | 59.30… | -4.10… | 1.03… | 130.81… | 1.74… | -1.75… | 0.038
… | … | … | … | … | … | … | … | …
… | … | … | … | … | … | … | … | …

## References

  * Ahn et al. (2012) Ahn, C. P., Alexandroff, R., Allende Prieto, C., et al. 2012, ApJS, 203, 21
  * Alfaro-Cuello et al. (2019) Alfaro-Cuello, M., Kacharov, N., Neumayer, N., et al. 2019, ApJ, 886, 57
  * Allende Prieto et al. (2006) Allende Prieto, C., Beers, T. C., Wilhelm, R., et al. 2006, The Astrophysical Journal, 636, 804.  https://dx.doi.org/10.1086/498131
  * Allende Prieto et al. (2020) Allende Prieto, C., Cooper, A. P., Dey, A., et al. 2020, Research Notes of the American Astronomical Society, 4, 188
  * Allende Prieto et al. (2023) Allende Prieto, C., Aguado, D. S., González Hernández, J. I., et al. 2023, ApJ, 957, 76
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167
  * Awad et al. (2024) Awad, P., Canducci, M., Balbinot, E., et al. 2024, A&A, 683, A14
  * Banik et al. (2018) Banik, N., Bertone, G., Bovy, J., & Bozorgnia, N. 2018, J. Cosmology Astropart. Phys, 2018, 061
  * Banik et al. (2021) Banik, N., Bovy, J., Bertone, G., Erkal, D., & de Boer, T. J. L. 2021, J. Cosmology Astropart. Phys, 2021, 043
  * Belokurov et al. (2018) Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611
  * Belokurov & Kravtsov (2024) Belokurov, V., & Kravtsov, A. 2024, MNRAS, 528, 3198
  * Bennett & Bovy (2019) Bennett, M., & Bovy, J. 2019, MNRAS, 482, 1417
  * Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition
  * Bonaca et al. (2019) Bonaca, A., Hogg, D. W., Price-Whelan, A. M., & Conroy, C. 2019, ApJ, 880, 38
  * Bonaca et al. (2020) Bonaca, A., Conroy, C., Hogg, D. W., et al. 2020, The Astrophysical Journal Letters, 892, L37.  https://dx.doi.org/10.3847/2041-8213/ab800c
  * Bonaca et al. (2021) Bonaca, A., Naidu, R. P., Conroy, C., et al. 2021, ApJ, 909, L26
  * Bowden et al. (2015) Bowden, A., Belokurov, V., & Evans, N. W. 2015, MNRAS, 449, 1391
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Carlberg (2012) Carlberg, R. G. 2012, ApJ, 748, 20
  * Carlberg (2013) —. 2013, ApJ, 775, 90
  * Carlberg (2018) —. 2018, ApJ, 861, 69
  * Carlberg & Agler (2023) Carlberg, R. G., & Agler, H. 2023, ApJ, 953, 99
  * Carpenter et al. (2017) Carpenter, B., Gelman, A., Hoffman, M. D., et al. 2017, Journal of Statistical Software, 76, 1
  * Chen & Gnedin (2022) Chen, Y., & Gnedin, O. Y. 2022, MNRAS, 514, 4736
  * Chen & Gnedin (2023) —. 2023, MNRAS, 522, 5638
  * Choksi & Gnedin (2019) Choksi, N., & Gnedin, O. Y. 2019, MNRAS, 486, 331
  * Choksi et al. (2018) Choksi, N., Gnedin, O. Y., & Li, H. 2018, MNRAS, 480, 2343
  * Cooper et al. (2023) Cooper, A. P., Koposov, S. E., Allende Prieto, C., et al. 2023, ApJ, 947, 37
  * Creasey et al. (2019) Creasey, P., Sales, L. V., Peng, E. W., & Sameie, O. 2019, MNRAS, 482, 219
  * de Boer et al. (2018) de Boer, T. J. L., Belokurov, V., Koposov, S. E., et al. 2018, MNRAS, 477, 1893
  * de Boer et al. (2020) de Boer, T. J. L., Erkal, D., & Gieles, M. 2020, MNRAS, 494, 5315
  * DESI Collaboration et al. (2016) DESI Collaboration, Aghamousa, A., Aguilar, J., et al. 2016, arXiv e-prints, arXiv:1611.00036
  * DESI Collaboration et al. (2022) DESI Collaboration, Abareshi, B., Aguilar, J., et al. 2022, AJ, 164, 207
  * DESI Collaboration et al. (2023a) DESI Collaboration, Adame, A. G., Aguilar, J., et al. 2023a, arXiv e-prints, arXiv:2306.06308
  * DESI Collaboration et al. (2023b) —. 2023b, arXiv e-prints, arXiv:2306.06307
  * Dey et al. (2019) Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168
  * Dillamore et al. (2022a) Dillamore, A. M., Belokurov, V., Evans, N. W., & Price-Whelan, A. M. 2022a, MNRAS, 516, 1685
  * Dillamore et al. (2022b) Dillamore, A. M., Belokurov, V., Font, A. S., & McCarthy, I. G. 2022b, MNRAS, 513, 1867
  * Dodge et al. (2023) Dodge, B. C., Slone, O., Lisanti, M., & Cohen, T. 2023, MNRAS, 518, 2870
  * Doke & Hattori (2022) Doke, Y., & Hattori, K. 2022, ApJ, 941, 129
  * Drimmel & Poggio (2018) Drimmel, R., & Poggio, E. 2018, Research Notes of the American Astronomical Society, 2, 210
  * Earp et al. (2017) Earp, S. W. F., Debattista, V. P., Macciò, A. V., & Cole, D. R. 2017, MNRAS, 469, 4095
  * Erkal & Belokurov (2015) Erkal, D., & Belokurov, V. 2015, MNRAS, 450, 1136
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gaia Collaboration (2022) Gaia Collaboration. 2022, VizieR Online Data Catalog, I/355
  * Gialluca et al. (2021) Gialluca, M. T., Naidu, R. P., & Bonaca, A. 2021, The Astrophysical Journal Letters, 911, L32.  https://doi.org/10.3847%2F2041-8213%2Fabf491
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration, Abuter, R., Amorim, A., et al. 2018, A&A, 615, L15
  * Grillmair & Dionatos (2006) Grillmair, C. J., & Dionatos, O. 2006, ApJ, 643, L17
  * Hahn et al. (2023) Hahn, C., Wilson, M. J., Ruiz-Macias, O., et al. 2023, AJ, 165, 253
  * Han et al. (2023) Han, J. J., Conroy, C., & Hernquist, L. 2023, Nature Astronomy, 7, 1481
  * Han et al. (2022) Han, J. J., Naidu, R. P., Conroy, C., et al. 2022, ApJ, 934, 14
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357\.  https://www.nature.com/articles/s41586-020-2649-2
  * Helmi (2004) Helmi, A. 2004, ApJ, 610, L97
  * Helmi et al. (2018) Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Nature, 563, 85
  * Hogg et al. (2010) Hogg, D. W., Bovy, J., & Lang, D. 2010, arXiv e-prints, arXiv:1008.4686
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90
  * Ibata et al. (2020) Ibata, R., Thomas, G., Famaey, B., et al. 2020, arXiv e-prints, arXiv:2002.01488
  * Ibata et al. (2023) Ibata, R., Malhan, K., Tenachi, W., et al. 2023, arXiv e-prints, arXiv:2311.17202
  * Ibata et al. (2002) Ibata, R. A., Lewis, G. F., Irwin, M. J., & Quinn, T. 2002, MNRAS, 332, 915
  * Johnston et al. (2005) Johnston, K. V., Law, D. R., & Majewski, S. R. 2005, ApJ, 619, 800
  * Johnston et al. (2002) Johnston, K. V., Spergel, D. N., & Haydn, C. 2002, ApJ, 570, 656
  * Kirby et al. (2013) Kirby, E. N., Cohen, J. G., Guhathakurta, P., et al. 2013, ApJ, 779, 102
  * Koposov et al. (2024) Koposov, S., Allende Prieto, C., Cooper, A. P., & DESI Collaboration. 2024, MNRAS submitted
  * Koposov (2019) Koposov, S. E. 2019, RVSpecFit: Radial velocity and stellar atmospheric parameter fitting, Astrophysics Source Code Library, record ascl:1907.013, Astrophysics Source Code Library, ascl:1907.013
  * Koposov et al. (2010) Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260
  * Koposov et al. (2019) Koposov, S. E., Belokurov, V., Li, T. S., et al. 2019, MNRAS, 485, 4726
  * Koposov et al. (2023) Koposov, S. E., Erkal, D., Li, T. S., et al. 2023, MNRAS, 521, 4936
  * Küpper et al. (2015) Küpper, A. H. W., Balbinot, E., Bonaca, A., et al. 2015, ApJ, 803, 80
  * Küpper et al. (2012) Küpper, A. H. W., Lane, R. R., & Heggie, D. C. 2012, MNRAS, 420, 2700
  * Law & Majewski (2010) Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229
  * Law et al. (2009) Law, D. R., Majewski, S. R., & Johnston, K. V. 2009, ApJ, 703, L67
  * Levi et al. (2013) Levi, M., Bebek, C., Beers, T., et al. 2013, arXiv e-prints, arXiv:1308.0847
  * Li et al. (2018) Li, G.-W., Yanny, B., & Wu, Y. 2018, ApJ, 869, 122
  * Li et al. (2022) Li, T. S., Ji, A. P., Pace, A. B., et al. 2022, ApJ, 928, 30
  * Lindegren et al. (2021a) Lindegren, L., Bastian, U., Biermann, M., et al. 2021a, A&A, 649, A4
  * Lindegren et al. (2021b) Lindegren, L., Klioner, S. A., Hernández, J., et al. 2021b, A&A, 649, A2
  * Malhan & Ibata (2019) Malhan, K., & Ibata, R. A. 2019, MNRAS, 486, 2995
  * Malhan et al. (2019a) Malhan, K., Ibata, R. A., Carlberg, R. G., et al. 2019a, ApJ, 886, L7
  * Malhan et al. (2019b) Malhan, K., Ibata, R. A., Carlberg, R. G., Valluri, M., & Freese, K. 2019b, ApJ, 881, 106
  * Malhan et al. (2018) Malhan, K., Ibata, R. A., Goldman, B., et al. 2018, MNRAS, 478, 3862
  * Malhan et al. (2020) Malhan, K., Valluri, M., & Freese, K. 2020, Monthly Notices of the Royal Astronomical Society, 501, 179.  https://doi.org/10.1093/mnras/staa3597
  * Malhan et al. (2022a) Malhan, K., Valluri, M., Freese, K., & Ibata, R. A. 2022a, ApJ, 941, L38
  * Malhan et al. (2022b) Malhan, K., Ibata, R. A., Sharma, S., et al. 2022b, ApJ, 926, 107
  * Marigo et al. (2017) Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77
  * Marino et al. (2011) Marino, A. F., Villanova, S., Milone, A. P., et al. 2011, ApJ, 730, L16
  * Meng & Gnedin (2022) Meng, X., & Gnedin, O. Y. 2022, MNRAS, 515, 1065
  * Miller et al. (2023) Miller, T. N., Doel, P., Gutierrez, G., et al. 2023, arXiv e-prints, arXiv:2306.06310
  * Myeong et al. (2019) Myeong, G. C., Vasiliev, E., Iorio, G., Evans, N. W., & Belokurov, V. 2019, MNRAS, 488, 1235
  * Ngan et al. (2015) Ngan, W., Bozek, B., Carlberg, R. G., et al. 2015, ApJ, 803, 75
  * Ngan et al. (2016) Ngan, W., Carlberg, R. G., Bozek, B., et al. 2016, ApJ, 818, 194
  * Nibauer et al. (2023) Nibauer, J., Bonaca, A., Lisanti, M., Erkal, D., & Hastings, Z. 2023, arXiv e-prints, arXiv:2312.09233
  * Nogueras-Lara (2022) Nogueras-Lara, F. 2022, A&A, 666, A72
  * Pastorelli et al. (2019) Pastorelli, G., Marigo, P., Girardi, L., et al. 2019, MNRAS, 485, 5666
  * Pastorelli et al. (2020) —. 2020, MNRAS, 498, 3283
  * Pedregosa et al. (2011) Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning Research, 12, 2825
  * Price-Whelan et al. (2022) Price-Whelan, A., Sipőcz, B., Wagg, T., et al. 2022, adrn/gala: v1.6.1, Zenodo, vv1.6.1, Zenodo, doi:10.5281/zenodo.7299506
  * Price-Whelan (2017) Price-Whelan, A. M. 2017, The Journal of Open Source Software, 2, doi:10.21105/joss.00388.  https://doi.org/10.21105%2Fjoss.00388
  * Price-Whelan & Bonaca (2018) Price-Whelan, A. M., & Bonaca, A. 2018, The Astrophysical Journal Letters, 863, L20.  https://dx.doi.org/10.3847/2041-8213/aad7b5
  * Reid & Brunthaler (2004) Reid, M. J., & Brunthaler, A. 2004, ApJ, 616, 872
  * Sanders et al. (2016) Sanders, J. L., Bovy, J., & Erkal, D. 2016, MNRAS, 457, 3817
  * Schlafly et al. (2010) Schlafly, E. F., Finkbeiner, D. P., Schlegel, D. J., et al. 2010, ApJ, 725, 1175
  * Siegal-Gaskins & Valluri (2008) Siegal-Gaskins, J. M., & Valluri, M. 2008, ApJ, 681, 40
  * Silber et al. (2023) Silber, J. H., Fagrelius, P., Fanning, K., et al. 2023, AJ, 165, 9
  * Simon (2019) Simon, J. D. 2019, ARA&A, 57, 375
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287
  * Vasiliev et al. (2021) Vasiliev, E., Belokurov, V., & Erkal, D. 2021, MNRAS, 501, 2279
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261\.  http://www.nature.com/articles/s41592-019-0686-2
  * Yan et al. (2022) Yan, H., Li, H., Wang, S., et al. 2022, The Innovation, 3, 100224
  * Zaritsky & Behroozi (2023) Zaritsky, D., & Behroozi, P. 2023, MNRAS, 519, 871