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

¹Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA; email:[email protected]

²NSF NOIRLab, 950 N. Cherry Ave., Tucson, AZ 85719, USA

³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

⁵Department of Astronomy & Astrophysics, University of Toronto, Toronto, ON M5S 3H4, Canada

⁶Department of Astronomy, University of Michigan, Ann Arbor, MI 48109, USA

⁷Institute of Astronomy and Department of Physics, National Tsing Hua University, 101 Kuang-Fu Rd. Sec. 2, Hsinchu 30013, Taiwan

⁸Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

⁹Physics Dept., Boston University, 590 Commonwealth Avenue, Boston, MA 02215, USA

¹⁰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

¹²Institute of Astronomy, Madingley Rd, Cambridge CB3 0HA, UK

¹³Steward Observatory, University of Arizona, 933 N, Cherry Ave, Tucson, AZ 85721, USA

¹⁴Department of Physics & Astronomy, University College London, Gower Street, London, WC1E 6BT, UK

¹⁵Department of Physics and Astronomy, The University of Utah, 115 South 1400 East, Salt Lake City, UT 84112, USA

¹⁶Departamento de Física, Universidad de los Andes, Cra. 1 No. 18A-10, Edificio Ip, CP 111711, Bogotá, Colombia

¹⁷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

²⁰Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

²¹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

²⁴Sorbonne Université, CNRS/IN2P3, Laboratoire de Physique Nucléaire et de Hautes Energies (LPNHE), FR-75005 Paris, France

²⁵Instituto de Física, Universidad Nacional Autónoma de México, Cd. de México C.P. 04510, México

²⁶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

²⁸Department of Astronomy, The Ohio State University, 4055 McPherson Laboratory, 140 W 18th Avenue, Columbus, OH 43210, USA

²⁹Institució Catalana de Recerca i Estudis Avançats, Passeig de Lluís Companys, 23, 08010 Barcelona, Spain

³⁰Department of Physics and Astronomy, Siena College, 515 Loudon Road, Loudonville, NY 12211, USA

³¹Department of Physics & Astronomy, University of Wyoming, 1000 E. University, Dept. 3905, Laramie, WY 82071, USA

³²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

³⁴Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía, s/n, E-18008 Granada, Spain

³⁵Department of Physics, Kansas State University, 116 Cardwell Hall, Manhattan, KS 66506, USA

³⁶Department of Physics and Astronomy, Sejong University, Seoul, 143-747, Korea

³⁷Institute for Computational Cosmology, Department of Physics, Durham University, South Road, Durham DH1 3LE, UK $^{38}$CIEMAT, Avenida Complutense 40, E-28040 Madrid, Spain

$^{39}$Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA

$^{40}$University of Michigan, Ann Arbor, MI 48109, USA

$^{41}$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

$^{42}$Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, Menlo Park, CA 94305, USA

$^{43}$Physics Department, Stanford University, Stanford, CA 93405, USA

$^{44}$SLAC National Accelerator Laboratory, Menlo Park, CA 94305, USA

$^{45}$National Astronomical Observatories, Chinese Academy of Sciences, A20 Datun Rd., Chaoyang District, Beijing, 100012, P.R. China

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 [56]; Johnston et al. [63]; Law et al. [73]; Law & Majewski [72]; Koposov et al. [66]; Malhan & Ibata [79]; Vasiliev et al. [109]). 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 ($\lesssim2\ \mathrm{km\ s^{-1}}$), unperturbed tidal streams produced by GCs are both thin and kinematically cold. Therefore velocity dispersion variations, density perturbations and gaps along a GC stream might hold clues to past flyby 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. [64]; Ibata et al. [62]; Siegal-Gaskins & Valluri [105]; Carlberg [21] (2013); Erkal & Belokurov [46]; Bonaca et al. [15]).

The GD-1 stellar stream, first identified by Grillmair & Dionatos [51], likely originated as a GC (Koposov et al. [66]). 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. [66]; Bowden et al. [18]; Malhan & Ibata [79]) 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. [33]; Malhan et al. [82]; Price-Whelan & Bonaca [101]). 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 [101]; de Boer et al. [34]) which have been claimed as evidence of possible impacts by dark matter subhalos (Bonaca et al. [15]). Since the stream does not pass through the inner disk region (it has a perigalactic (apogalactic) distance of $\sim14$ (26) kpc, and inclination of 39° to the disk plane, Koposov et al. [66]) 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 [43]).

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. [71, 70]). 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. [60] conclude that the density profile of the stream shows periodic gaps separated by $\sim2.64 \pm0.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. [81] used Gaia DR2 data to show that the GD-1 stellar stream possesses a secondary diffuse and extended stellar component ($\sim100$ 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; Malhan et al. 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 [24]) and controlled simulations (Malhan et al. [83]). 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; Malhan et al. [84]), 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 ($\gtrsim4 \times10^{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 $\sim1 - 2$ km s$^{-1}$ ([16, 49, 60]). An additional $\sim200$ spectroscopically confirmed members are available from the SDSS and LAMOST surveys but have lower radial velocity precision ($\sim5 - 20$ km s$^{-1}$). In this work we present 115 new spectroscopically confirmed members of the GD-1 stream obtained with the Dark Energy Spectroscopic Instrument (DESI) with $\sim1 - 10$ km s$^{-1}$ velocity precision (median precision $2.7$ km s$^{-1}$) and $\sim0.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 ([37, 69]). Additional data releases are planned at regular intervals. In § 1.1, we briefly describe the DESI Survey Validation data set in the region overlapping GD-1. In § 2 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 § 3 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 § 4 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.

DESI OBSERVATIONS OF THE GD-1 REGION

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 $\sim3^\circ$ field of view ([35, 106, 90]). 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 ([35]).

DESI was installed in 2019 and commissioning and survey validation extended through May 2021 ([37]). 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) ([35, 74]).

Over the first 5-years of the DESI Milky Way Survey ([4, 31]) it is expected that spectra of $\simeq7$ 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. [39]). 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. [31]). DESI-MWS is expected to achieve an average spectroscopic completeness of $\sim28\%$ 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. [52]), 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 publicly 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.

SV3. The MWS target selection for SV was previously described in Allende Prieto et al. [4]. Other DESI SV programs are discussed in DESI Collaboration et al. [37]. 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. [31], 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. [69]).

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 a coordinate transformation described in Koposov et al. [66]¹, which also yields proper motion in $[\phi_1,\phi_2]$, hereafter referred to as $\mu_{\phi_1}$ and $\mu_{\phi_2}$² respectively. Most targets were also identified in the DESI Legacy Survey (DECaLS) (Dey et al. [39]), from which the $g$ and $r$-band magnitudes for these stars were measured. The Legacy Survey magnitudes were dereddened using dust maps [104] 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 [48]. 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’) [102, 44, 13] (GRAVITY Collaboration et al. 2018). Distances to individual stars were estimated using the distance modulus in Equation (1) below (also see Appendix § 5.1). Gaia parallaxes were zero-point corrected following the prescriptions in [77].

The 2D histogram in Figure 1 shows the surface density of Gaia DR3 stars with DECaLS $0.2 \le(g-r) < 0.4$ and $18 \le r < 20$, $-10\,\mathrm{mas}\,\mathrm{yr}^{-1} < \mu_{\phi_1} < -4\,\mathrm{mas}\,\mathrm{yr}^{-1}$ and $|\mu_{\phi_2}| < -5\,\mathrm{mas}\,\mathrm{yr}^{-1}$. These selections are similar to and based on previous works (e.g., [101]). 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$.

The stellar spectra obtained by the DESI spectrograph are processed with the RVS pipeline, based on RVSpec [65], to measure each star’s line-of-sight (radial) velocity and stellar atmospheric parameters ($T_{\mathrm{eff}}$, $\log(g)$, [Fe/H], [$\alpha$/Fe]). Using the radial velocities from RVS, the spectra are independently processed by the SP pipeline which is based on FERRE [3].³ In addition to the 4 atmospheric parameters, above the SP pipeline measures the abundances of Mg, C, Ca.

The RVS and SP pipeline outputs have been calibrated against the APOGEE measurements. The DESI EDR dataset contains spectra for 507 stars that are also in the APOGEE DR17 catalog [69]. 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 $[\mathrm{Fe}/\mathrm{H}]_{\mathrm{DESI}}-[\mathrm{Fe}/\mathrm{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 pipeline. For further details of the pipelines and their validation see Cooper et al. [31] and Koposov et al. [69]. 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] in this work.

Despite the good performance of both RVS and SP pipelines on the 507 APOGEE stars however, Koposov et al. [69] point out that the SP pipeline output is extremely sensitive to the S/N of the spectra, $T_{\mathrm{eff}}$, and discontinuities in the stellar templates, and therefore several quality cuts are recommended to obtain robust results (see Appendix § 5.1 of Koposov et al. [69] for details). We find that less than 50% of the stars in the GD-1 sample pass the quality cuts needed for robust SP abundances but all pass the RVS cuts. Therefore for the rest of this paper we restrict ourselves to using the output of the RVS pipeline.

GD-1 Stream Membership

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 [101] and [79].

The distance to the GD-1 stream is known to vary with $\phi_1$ [66]. We assume a distance modulus $\mathcal{D}$ given by,

which is similar to the function estimated by [66]. This relation is a simple quadratic fit to the subgiant branch stars in GD-1. More details are provided in Appendix § 5.1 where we also compare this relation with the distance relations proposed by several other authors [101, 33] (Li et al. 2018). The functions defined by de Boer et al. [33] and Li et al. (2018) differ only slightly from the distance-$\phi_1$ relation in Equation (1) (see Figure 13 in Appendix § 5.1). We use equation 1 in the rest of the paper.

We impose a Gaia parallax cut such that

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_{\mathrm{err}}$ is the uncertainty on the Gaia parallax for each star, and $\varpi_{\mathrm{sys}} = 0.05\,\mathrm{mas}$ is the systematic error on the parallax [78].

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). With the above selections we compute $\delta\mu_{\phi_1}$ and $\delta\mu_{\phi_2}$ relative to the proper motion stream tracks. The densest part of the 2D distribution in $\delta\mu_{\phi_1}$ and $\delta\mu_{\phi_2}$ is reasonably well fitted by a 2D circular Gaussian of standard deviation 0.5 mas/yr. We conservatively select stars within 5 times this standard deviation i.e. all stars within $\pm2$ mas/yr around each proper motion spline to obtain the selection shown by black points in Figure 2 (c) which is similar to the polygon defined in [101]. 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 § 2.3) are given in Table 1.

$\phi_1$	$\mu_{\phi_1}$	$\mu_{\phi_2}$	$\phi_2$	$V_{\mathrm{GSR}}$
[°]	[mas/yr]	[mas/yr]	[°]	[km/s]
-90.0	-3.00	-4.32	-0.70	75.49
-70.0	-1.50	-6.92	-0.88	60.78
-60.0	-0.72	-7.46	-0.70	29.15
-46.0	-0.10	-7.54	-0.37	-21.31
-37.6	0.05	-7.29	-0.21	-46.37
-29.2	0.22	-6.85	-0.09	-66.46
-20.8	0.21	-6.28	0.00	-93.05
-12.4	-0.18	-5.61	0.08	-113.76
-4.0	-0.84	-4.90	0.16	-133.84
10.0	-2.03	-3.70	0.32	-189.77

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. [20]; Tang et al. [108]; Marigo et al. [86]; Pastorelli et al. [96, 97]).⁴ Other authors have used similar isochrones selections (e.g., [60]. We accepted stars along the PARSEC isochrone if they lie within $\delta(g-r) \pm0.1$ (shown by faint red curves).

Radial Velocity Distribution

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_{\mathrm{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| \leq3^\circ$. This value of $|\delta\phi_2|$ is motivated by the fact that in EDR, DESI only observed within $3^\circ$ from the stream track. Future data releases will include much wider coverage. Stars observed by DESI from the EDR catalog

[37] are shown as black points. Stars observed by SDSS ([1]) (blue) DR9⁵ and LAMOST DR8⁶ are shown by blue symbols. The 43 stars identified as GD-1 members using MMT-Hectochelle spectra ([16]) are shown by orange symbols. The selection of ([16]) only considered stars with velocities within $7$ km s$^{-1}$ from the stream track or the spur as GD-1 members, but we do not expand on this selection since few 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 § 5.2.

TARGET-ID	RA	Dec	$\phi_1$	$\phi_2$	$V_{\mathrm{GSR}}$	$V_{\mathrm{err}}$	[Fe/H]	$[\mathrm{Fe/H}]_{\mathrm{err}}$
[deg]	[deg]	[deg]	[deg]	[deg]	[km s$^{-1}$]	[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.04
39633377...	191.69...	59.30...	-4.10...	1.03...	130.81...	1.74...	-1.75...	0.04
...	...	...	...	...	...	...	...	...
...	...	...	...	...	...	...	...	...

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 [16] 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 115 new spectroscopically confirmed radial velocities and improved RV measurements for 17 previously observed members.

A spline fit to the stream track is shown with a red curve (partially solid, partially dotted) (see below for details of how this fit is obtained) along with dashed lines showing $\delta V_{\mathrm{GSR}} = \pm30$ km s$^{-1}$ from the stream track. Figure 4 (b) shows kernel density estimates of $\delta V_{\mathrm{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. 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_{\mathrm{GSR}}| \leq7$ km s$^{-1}$ imposed by [16].

In order to obtain the stream track, we use a mixture modeling approach that properly accounts for the uncertainties in the radial velocity measurements and background contamination. 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 \leq0^\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 $v$ conditional on $\phi_1$ is given by,

which is identical to equation 2 of Koposov et al. [68]. In the above equation $f(\phi_1)$ is the mixing fraction, $V(\phi_1)$ is the velocity of the stream, and $S_V(\phi_1)$ is the velocity dispersion of the stream, all of which are represented as cubic splines with the values at the knots being model parameters. $v_{\mathrm{bg}}$ is the mean velocity of the Gaussian describing the background, $d_{\mathrm{bg}}$ is the gradient of the mean velocity of the background with respect to $\phi_1$, and $S_{\mathrm{bg}}$ is the velocity dispersion of the background.

The model is implemented using the Stan probabilistic programming language [25] and is sampled using the 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. [68]. 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. The 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. We note that although the method yields a measurement of the velocity dispersion along the stream, the main purpose of this analysis is to robustly obtain the stream track in the presence of the background.

Figure 5 (top) shows the measured radial velocity $V_{\mathrm{GSR}}$ vs. $\phi_1$ for all DESI targets within $\pm3^\circ$ from the stream track (grey points). The uncertainties on the velocities are not visible as they are $\lesssim15$ km/s. Red lines show posterior samples from radial velocity curves. Spline knots (obtained with Stan) are shown as black dots.

The lower panel of Figure 5 shows the velocity dispersion along the stream in the region $-46^\circ\leq\phi_1 \leq-4^\circ$. The error-bars show the $1\sigma$ uncertainties defined as the 16th-84th percentile credible intervals determined from the posterior samples (also given in Table 2). The velocity dispersion everywhere is between 2.36 and 5.47 km s$^{-1}$ with the uncertainties of $\pm0.7$–$\pm2.5$ km s$^{-1}$.

$\phi_1$	$\sigma_{V_{\mathrm{GSR}}}$	$1\sigma$ uncertainty
degrees	km s$^{-1}$	km s$^{-1}$
-46.0	2.36	1.24
-32.0	5.47	1.26
-18.0	3.78	0.73
-4.0	2.72	2.45

We will return to the issue of the radial velocity dispersion in § 3.

The Stan spline fit above is obtained by using only DESI velocity measurements (since these have the smallest uncertainties) and is shown by the solid red curve in Figure 4. Beyond the range of $\phi_1$ values modeled by the Stan spline, we select knots by eye, rather than by using the median/mean velocity in bins primarily because there are few DESI detections in these regions and the outliers are not symmetrically distributed relative to the stream track (see Figure 5). The 3-year DESI sample has substantially more observations available and will allow for a more robust analysis along the entire length of the stream.

Metallicity distribution of stream stars

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 § 2.1 and § 2.2. Stars with the additional constraint that they lie within 30 km s$^{-1}$ from the stream track in Figure 4 (§ 2.3) are shown in black. It is clear from the black histogram in Figure 6 that only one star has [Fe/H] $> -1.5$ and one has [Fe/H] $< -3.5$. We exclude these stars bringing the sample from DESI to 151 (115 new, 36 improved). With these metallicity cuts 186 GD-1 stars from SDSS+LAMOST remain, along with 43 stars from MMT ([16]).

!Histogram of the [Fe/H] distribution for stars in the GD-1 stream bbox=0.13,0.11,0.45,0.43

As discussed in § 1.1 comparisons with 507 DESI stars in also observed by APOGEE indicate that the RVS metallicity has a median offset [Fe/H]$_{\mathrm{DESI}}$−[Fe/H]$_{\mathrm{APOGEE}} = 0.06^{+0.13}_{-0.18}$ dex. However, despite this small offset, we note that Cooper et al. [31] showed that [Fe/H] measurements obtained with DESI for several GCs (M13, M92, NGC 5053) showed significant spreads of 0.3 dex or greater⁷. Furthermore they often had mean metallicities that were systematically lower by 0.2 to 0.4 dex relative to literature values. Given these rather large observed spread of 0.2-0.3dex in the [Fe/H] values recovered by the RVS pipeline we caution readers against inferring too much from the spread in $[\mathrm{Fe}/\mathrm{H}]$ for GD-1 members reported in this paper.

Despite these concerns with current DESI metallicity measurements the mean $[\mathrm{Fe}/\mathrm{H}] = -2.49 \pm0.25$ for all stars in the black histogram is consistent with values from previous works: $[\mathrm{Fe}/\mathrm{H}] = -2.24 \pm0.21$ ([79]), $[\mathrm{Fe}/\mathrm{H}] = -2.3 \pm0.1$ ([16]), and $[\mathrm{Fe}/\mathrm{H}] = -2.49 \pm0.03$ Martin et al. [88] and $[\mathrm{Fe}/\mathrm{H}] = -2.47$ ([61]).

Identification and Characterization of the Cocoon

As discussed in § 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 § 4), 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| \le3^\circ$), can be seen in Figure 7.

Figure 7. This plot only shows observations with full 6-D phase space and metallicity information. Malhan et al. [81] identified the GD-1 cocoon primarily in the region $-10^\circ< \phi_1 < 0^\circ$ based 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 $-30^\circ< \phi_1 < 2^\circ$.

We now analyze the phase space distribution of the stars in the region bounded 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 to larger distances perpendicular to the stream are required to confidently assess the 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 ([101, 15]).

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 ($\pm30\ \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_{\mathrm{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_{\mathrm{GSR}}$ and scatter plots of $\delta\mu_{\phi_1}$ and $\delta\mu_{\phi_2}$ vs. $\delta V_{\mathrm{GSR}}$. We also note that the plots show evidence for a greater spread in $\delta V_{\mathrm{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). In particular we note that the significant spread in $V_{\mathrm{GSR}}$ at small $|\delta\mu_{\phi_1}|$, likely arises from cocoon stars in front of and behind the thin stream, evidence that the cocoon is probably a 3D cylindrical structure around the thin stream.

Multi-Gaussian Analysis

We wish to separate the “thin stream” and “cocoon” components in order to individually 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$.

fraction $Q$. We adapt the formalism of Hogg et al. [58] 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”). Assuming each component is well described by a Gaussian, the likelihood function can be written as,

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 (5) 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. [47]). 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_{\mathrm{corr}} \approx80$ steps, we discard the first 1600 samples as “burn-in” and thin the remaining chain by a factor $N_{\mathrm{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 for 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

Parameter	Univariate		Bivariate
	Thin	Cocoon	Thin	Cocoon
$\overline{\delta\phi_2}$ [deg]	$-0.118 \pm0.017$	$-0.149 \pm0.161$	$-0.111 \pm0.017$	$-0.161 \pm0.157$
$\Sigma_{\delta\phi_2}$ [deg]	$0.149 \pm0.014$	$1.265 \pm0.122$	$0.149 \pm0.014$	$1.252 \pm0.121$
$\overline{\delta V_{\mathrm{GSR}}}$ [km s$^{-1}$]	$\ldots$	$\ldots$	$-1.19 \pm0.52$	$1.67 \pm1.21$
$\Sigma_{\delta V_{\mathrm{GSR}}}$ [km s$^{-1}$]	$\ldots$	$\ldots$	$3.090 \pm0.76$	$7.81 \pm1.01$
$Q_{\mathrm{thin}}$	$0.63 \pm0.05$	$\ldots$	$0.61 \pm0.05$	$\ldots$

parameters in this table are $1\sigma$ uncertainties determined from the 16th and 84th percentile credible intervals of the posterior probability distributions.

Next we jointly fit the $\phi_2$ and $V_{\mathrm{GSR}}$ distributions with two bivariate-Gaussians described by the likelihood function,

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_{\mathrm{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,\mathrm{thin}} = (0.149 \pm0.014)^\circ$ implies a FWHM of the “thin stream” component in $\phi_2$ to be $2.355 \times\sigma_{\phi_2,\mathrm{thin}} = 0.351^\circ\approx55$ pc (at a mean distance of 9 kpc) while the cocoon component is much broader with standard deviation $\Sigma_{\delta\phi_2,\mathrm{cocoon}} = (1.252 \pm0.121)^\circ$

or FWHM $\approx2.948^\circ\approx463$ pc. For the distribution in $V_{\mathrm{GSR}}$ the bivariate analysis gives velocity dispersions $\Sigma_{V_{\mathrm{GSR},\mathrm{thin}}} = (3.090 \pm0.76)$ km s$^{-1}$ and $\Sigma_{V_{\mathrm{GSR},\mathrm{cocoon}} = (7.81 \pm1.01)}$ km s$^{-1}$ for the thin stream and cocoon components respectively, with the fraction of stars in the thin stream $Q = 0.61 \pm0.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_{\mathrm{GSR}}$ yields nearly identical results for the angular widths of the thin stream and cocoon. Both methods also give $Q \sim0.62$ 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_{\mathrm{thin}} \ge0.5$ as determined by the bivariate analysis considering both position and velocity $\delta V_{\mathrm{GSR}}$. The vertical lines are located at $\overline{\delta\phi}_{2,\mathrm{thin}} \pm\Sigma_{\delta\phi_{2,\mathrm{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}$. As stated previously the spread in $\delta V_{GSR}$ could be due to high velocity cocoon stars in front of and behind the thin component of the stream. In addition, over this range of $\phi_1$, the distance gradient implies that the stream is at an angle of $\sim30^\circ$ to our line-of-sight. Thus the increased dispersion in $\delta V_{GSR}$ and $\delta\mu_{\phi_1}$ but not in $\delta\mu_{\phi_2}$ suggests that there has been greater heating along the stream than perpendicular to the stream.

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 § 3.1 to define each component. We then model each component’s velocity distribution with a univariate Gaussian for that component and represent the background by a much broader Gaussian in velocity truncated at $\delta V_{GSR} = \pm30$ km s$^{-1}$ (the velocity cut shown by the dashed lines in Figure 4a). The mean velocity and dispersion of each component (thin stream or cocoon) and the background contamination are fitted parameters. The likelihood function is similar to that given in Equation (5). Therefore we repeat the exercise of modeling a main component and a broad truncated Gaussian contaminant twice, once for the “thin stream” members and again for the “cocoon” members.

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

When this analysis is applied to the members identified with the physically “thin” component (stars colored black in Figure 10 the model yields a thin stream velocity dispersion of $2.85 \pm0.6$ km s$^{-1}$ and the fraction of stars assigned to the background was $0.12 \pm0.05$. This velocity dispersion for the thin stream component is consistent (within the uncertainties) to the dispersion obtained from the bivariate multi-Gaussian analysis in the previous section.

When applied to the members of the physically broad cocoon component (stars colored red in Figure 10 we find that the cocoon has a smaller velocity dispersion of $4.84 \pm1.13$ km s$^{-1}$ and the background fraction in the cocoon is also 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 Figure 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 $\sim2.85 - 3.09$ km s$^{-1}$, while the cocoon has a velocity dispersion between $4.84$-$7.81$ 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.

We remind the reader that the Stan analysis in section 3.3, which was primarily intended to derive the stream track in $V_{GSR}-\phi_1$ space while properly accounting for background contamination, also resulted in velocity dispersions in the $-32 \leq\phi_1 \leq-4$ of $5.47 \pm1.26$ to $2.72 \pm2.45$ km s$^{-1}$, which is consistent with the above estimates (see Table 2).

Our measurement of the dispersion of the ‘thin’ stream component ($\sigma_{V_{GSR,\mathrm{thin}}} \sim3$ km s$^{-1}$) is consistent with previous estimates of both other GC streams and the GD-1 stream itself ((Galiullina et al. 2021)). In contrast, the larger velocity dispersion of the broader cocoon component of $\sigma_{V_{GSR,\mathrm{cocoon}}} \sim5 - 8$ km s$^{-1}$ is comparable to the velocity dispersion of classical dwarf spheroidal galaxies like Carina, Draco, Fornax, Leo I etc. ([107]). 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$ km s$^{-1}$ with a mean value of around $11$ km s$^{-1}$ and are significantly larger than the dispersions of known GC streams (e.g., (Li et al. 2022)).

Metallicity of thin stream and cocoon components

Another way to assess if the cocoon velocity dispersion in § 3.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 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 and the metallicity of the cocoon stars is indistinguishable from that of the thin component within the uncertainties.

!Kernel density estimates of [Fe/H] for the thin stream and cocoon components bbox=0.17,0.10,0.45,0.29

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 $\langle[\mathrm{Fe}/\mathrm{H}]_{\mathrm{thin}}\rangle= -2.54 \pm0.04$ dex and intrinsic metallicity dispersion is $\sigma_{[\mathrm{Fe}/\mathrm{H}],\mathrm{thin}} = 0.27 \pm0.03$ dex.

The cocoon component has a mean $\langle[\mathrm{Fe}/\mathrm{H}]_{\mathrm{cocoon}}\rangle= -2.47 \pm0.06$ dex with dispersion $\sigma_{[\mathrm{Fe}/\mathrm{H}],\mathrm{cocoon}} = 0.31 \pm0.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. [61] 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.255$, 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 intrinsic metallicity dispersion of GD-1 is significantly higher than typical GC streams which have metallicity dispersions of $0.03$-$0.05$ dex. This rather high inferred intrinsic metallicity dispersion of GD-1 may be explained by the current systematic floor of $0.3$ dex in the uncertainty of [Fe/H] measurements resulting from the RV pipeline, as discussed in § 2.4 and as shown in Cooper et al. [31]. [Fe/H] measurements in GCs obtained with DESI show significantly larger spreads of $0.3$ dex or greater 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).

Discussion and conclusions

Summary of results

We report spectroscopic analysis of the GD-1 tidal stream using observations with Dark Energy Spectroscopic Instrument (DESI Collaboration et al. [35] (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 $[\mathrm{M}/\mathrm{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_{\mathrm{GSR}}$ coordinates using a Gaussian mixture modelling approach, accounting for contamination by a background population, following the approach described in Koposov et al. [68]. Rejecting one star 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 151 of which 115 are new members and for the remaining 36 stars, 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 $-30^\circ\leq\phi_1 \leq2^\circ$ of the GD-1 tidal stream (where the majority of our observations were obtained). This cocoon was first identified by Malhan et al. [81] using primarily Gaia positions and proper motions and was recently confirmed by Ibata et al. [61].

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 \pm0.014^\circ$ and $\Sigma_{\delta\phi_2,\mathrm{cocoon}} = 1.252 \pm0.121^\circ$. Since our DESI fields are restricted to $\lvert\delta\phi_2\rvert\sim\pm3^\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_{\mathrm{GSR},\mathrm{thin}}} = 3.09 \pm0.76\ \mathrm{km\ s^{-1}}$ and $\Sigma_{V_{\mathrm{GSR},\mathrm{cocoon}}} = 7.81 \pm1.01\ \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 $\sim4.84\ \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 $\langle[\mathrm{Fe}/\mathrm{H}]_{\mathrm{thin}}\rangle= -2.54 \pm0.04$ dex and $\langle[\mathrm{Fe}/\mathrm{H}]_{\mathrm{cocoon}}\rangle= -2.47 \pm0.06$ dex respectively, with intrinsic metallicity dispersions $\sigma_{[\mathrm{Fe}/\mathrm{H}],\mathrm{thin}} = 0.27 \pm0.03$ dex and $\sigma_{[\mathrm{Fe}/\mathrm{H}],\mathrm{cocoon}} = 0.31 \pm0.04$ 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.

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. 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. 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. A thin tidal stream generated from the GC progenitor could have been heated by multiple interactions with dark subhalos in the MW.

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 [89]. 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; Malhan et al. 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. [83] 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. [83] 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 a tidal stream with a velocity dispersion of $\sim2$–$5\ \mathrm{km\ s^{-1}}$, while a dwarf galaxy with a cuspy halo would produce a stream velocity dispersion of order 8-12 km $s^{-1}$. In this pre-accretion tidal stripping scenario, a cocoon velocity dispersion of $\sim8$ km $s^{-1}$ would imply that GD-1’s parent dwarf galaxy had a cuspy dark matter halo while $\sim5$ 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 its metallicity are reasonably consistent with those of other Local Group dwarf galaxies. 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., [e.g. 2, 95]), 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 its 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’ ([91, 17, 85]). Accreted GCs which are about 60% of the GC population ([30, 29, 32]) tend to be on more radial orbits and tend to be found at larger radii ([27, 12] (Chen & Gnedin 2023)). 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 ([14]). 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 GCs with masses at or above $\sim2 \times10^4 M_\odot$ ([26, 87, 19]) but not in dwarf galaxies. Various estimates of the mass of the GD-1 progenitor exist in the literature: $2 \times10^4 M_\odot$ ([66]), $5 \times10^4 M_\odot$ ([49]), $7 \times10^4 M_\odot$ ([15]) and $1 \times10^5$ ([101]). While these values are at the lower limit of the mass range of GCs showing multiple stellar populations and a Na-O anti-correlation, one reason to be optimistic that future studies may show such an anti-correlation is that a recent study of high resolution spectra of GD-1 stars (including some in the spur and blob) has already detected variations in C and Mg ([8]), suggesting the presence of multiple populations in this system.

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., [64, 62, 105, 21, 46, 103, 9, 15, 10, 24] (Carlberg 2013)). Simulations provide insights into the origin 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 ([92, 93]) 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 ([62, 105]). 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 ([49]), 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. [34] 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 \times10^{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” ([101]) (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. [80].

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. [61] 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 ($\sim10\ \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 \pm1.1\ \mathrm{km\ s^{-1}}$ while the cocoon has a velocity dispersion $\sigma= 29.1 \pm6.1\ \mathrm{km\ s^{-1}}$ with a fraction of $0.73 \pm0.09$ in the thin component.

Our analysis in § 3 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. [61]. 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 \pm0.38\ \mathrm{km\ s^{-1}}$ and $7.68 \pm0.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. [61]. 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 ([83, 24]).

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.

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 ([11, 57]) could have resulted in the disk of the Milky Way becoming misaligned with the resulting triaxial dark matter halo ([41, 54, 53]) 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 $\sim10 - 20\ \mathrm{km\ s^{-1}\ kpc^{-1}}$ (e.g., [45, 42]). Recently Nibauer et al. [94] 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 $\sim20\ \mathrm{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\ \mathrm{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.

GD-1 data expected from future DESI data releases

Future DESI data releases will include many more GD-1 members. The first 3 years of the survey have already yielded over 600 spectroscopically confirmed members 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 $\sim15$ km s$^{-1}$). Current plans for the future include observing up to ninety $3^\circ$ 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 $\sim3000$ GD-1 stream members, making it one of the most thoroughly surveyed GC streams.

MV gratefully acknowledges support from NASA-ATP grants 80NSSC20K0509 and 80NSSC24K0938. 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 NOIR-Lab; 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 NOIR-Lab 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 Gaia,⁸ processed by the Gaia Data Processing and Analysis Consortium (DPAC⁹). 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. The University of Michigan is located on the traditional territory of the Anishinaabe people. We acknowledge, their contemporary and ancestral ties to the land and their contributions to the university.

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

Software: astropy (Astropy Collaboration et al. 2018, 2022), emcee (Foreman-Mackey et al. [47]), gala (Price-Whelan [100]; Price-Whelan et al. [99]), matplotlib (Hunter [59]), numpy (Harris et al. [55]), scipy (Virtanen et al. [110]), scikit-learn (Pedregosa et al. [98]), stan (Carpenter et al. [25]).

Facilities: KPNO: Mayall (DESI), Gaia, Blanco (DECam).

Data Availability: Data behind all figures in this paper are available as FITS files at: https://zenodo.org/records/11638329 DOI: 10.5281/zenodo.11638329.

Appendix

Measuring the distance stream track

To calibrate the distance gradient along the stream, we use the proper motion and parallax selection described in § 2 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)$. Figure 13 (left) 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 (right) shows various relations for the distance to the GD-1 stream as a function of $\phi_1$ from Price-Whelan & Bonaca [101] (blue dot-dashed), de Boer et al. [33] (green dashed), Li et al. (2018) (black solid) and equation 1 (this work, red solid).

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 uncertainties, $V_{err}$ and [Fe/H]$_{err}$ are from the RVS pipeline. The $V_{GSR}$ measurements have been corrected by $0.93\,\mathrm{km}\,\mathrm{s}^{-1}$, the systematic error on RV measurements reported in [69]. No corrections have been applied to [Fe/H] measurements. TARGET_ID is the DESI target ID, RA and Dec values are from Gaia-DR3.