Typeset using

Keywords: Galaxy, halo -Galaxy, stellar content -surveys -galaxies, formation -Galaxy, structure

The density variations in thin stellar streams may encode important information on the nature of the dark matter. For instance, if dark matter aggregates into massive sub-halos, these perturbers are expected to scatter stars out of dynamically cold stellar streams, possibly leading to detectable gaps in those structures. Here we re-examine the density variations in the GD-1 stream, using Gaia DR2 astrometry, Pan-STARRS photometry, together with high precision radial velocities measured with the CFHT/ESPaDOnS and VLT/UVES instruments and complemented with public radial velocity catalogs. We show that after correcting for projection effects, the density profile exhibits high contrast periodic peaks, separated by 2.64 ± 0.18 kpc. An N-body simulation is presented that reproduces this striking morphology with simple epicyclic motion in a smooth Galactic potential. We also discuss the reliability of measuring density variations using ground-based photometric surveys, and for the particular case of GD-1 we highlight some of the artifacts present in the Gaia DR2 catalog along its track. Massive dark subhalos do not appear to be required to explain the density clumping along GD-1.

INTRODUCTION

Corresponding author: Rodrigo Ibata rodrigo.ibata@astro.unistra.fr

A key prediction of Λ Cold Dark Matter (ΛCDM) cosmology is that galaxies reside within dark matter halos that are composed of a hierarchy of smaller "sub-halos". This hierarchy is expected to continue down in mass all the way to a limit set by the (as yet unknown) thermal free-streaming length of the dark matter particle (Springel et al. 2008). Many thousands of these subhalos are expected to orbit within our Galaxy, but only the most massive would contain some baryonic component that could render them directly observable (as satellite galaxies). So detecting the huge predicted population of completely dark sub-halos requires identifying their gravitational influence on photons or on observable baryonic structures.

One promising avenue to detect the sub-halo population is to analyse the morphology and flux-ratios of strongly-lensed quasar images. At present the evidence appears consistent with ΛCDM (Ritondale et al. 2019;Hsueh et al. 2019), but at low statistical significance.

It is interesting therefore to consider how the subhalos may influence stellar sub-structures of our Milky Way or of other nearby galaxies. Although the expected fully-dark sub-halos could be very massive (up to ∼ 10 8 M ), their large physical scale makes these bodies very "fluffy", and interactions with the baryonic components of a galaxy will be subtle. One therefore needs to identify some dynamical probes that respond in a measurable way to small perturbations of the acceleration field. This realisation led several groups to propose that the fragility of dynamically-cold star streams could be used as a means to explore the sub-clustering of the dark matter on sub-galactic scales (Ibata et al. 2002;Johnston et al. 2002;Mayer et al. 2002).

Heating from a sub-halo fly-by will increase the velocity dispersion in a stream, and given that these initially can be very cold (e.g., the one-dimensional velocity dispersion in the GD-1 stream is ∼ 1 km s -1 , Malhan & Ibata 2019), the influence of the sub-halo flyby may be detectable, in principle. However, the practical difficulty in realising such a measurement is that streams generally possess a very low density of stars that are bright enough to be measured with good precision, which makes the dynamical heating effect challenging to detect.

A promising alternative to measuring velocity dispersion variations (which would require obtaining highprecision line-of-sight kinematics to hundreds or thousands of stars in a stream) is instead to make use of the stream's spatial morphology. Carlberg (2012) showed that characteristic underdensities or "gaps" are formed after a close flyby of a massive perturber. Indeed, for the specific case of the GD-1 stream, Carlberg (2016) proposed that sub-halos could be responsible for the gaps on scales of ∼ 10 • that were detected in Sloan Digital Sky Survey (SDSS) maps of the system (Carlberg & Grillmair 2013).

Recently, de Boer et al. (2019) have remeasured the morphology of the GD-1 stream using the excellent astrometric data from the Second Data Release (DR2) of the Gaia mission (Lindegren et al. 2018;Gaia Collaboration et al. 2018) combined with Pan-STARRS photometry (Chambers et al. 2016). Banik et al. (2019) use these data to detect a power spectrum of density variations along the stream that they claim requires the presence of a population of perturbing sub-halos of mass 10 7 M to 10 9 M with a density that is within the uncertainties of ΛCDM predictions.

The present work aims to examine these very interesting claims, providing additional data and analysis of the GD-1 system. The layout of this paper is as follows. Section 2 presents an overview of the GD-1 system, whose properties we re-derive in Section 3, based on a clean sample of stars, including new radial velocity measurements. With these new constraints we present the density profile along the stream in Section 4, finding that the profile is substantially more peaked than found by Banik et al. (2019), and displays periodic overdensities. In Section 5 we change tack to attempt to quantify the reliability of ground-based photometric surveys in order to estimate how confidently the surface density of a highly contaminated structure can be measured. In Section 6 we present some simple models to interpret the observed density profile. Our simulations show that the density spikes can be modelled by the escape of stars at low velocity from a globular cluster that has now completely dissolved. Finally, with all these caveats in mind, we measure the power spectrum of the GD-1 stream in Section 7. Our conclusions are laid out in Section 8.

THE GD-1 STREAM

The GD-1 stream was discovered by Grillmair & Dionatos (2006) in the SDSS, where it appeared as a 63 • long structure in matched filter maps designed to reveal metal-poor populations similar to that of the globular cluster M13 ([Fe/H] = -1.53, Harris 2010). The stream lies in the North Galactic cap region in the direction away from the Galactic center. Follow-up medium-resolution spectroscopy obtained by Koposov et al. (2010) showed that GD-1 has a relatively circular but retrograde orbit, with a pericenter at 14 kpc and an apocenter at 26 kpc. This orbit keeps the system well away from the inner regions of the Galactic disk, where interactions with giant molecular clouds could cause additional heating (Amorisco et al. 2016), that could contaminate the sought-for signal from ΛCDM substructure.

The advent of the Gaia DR2 catalog enabled the search for streams over the full sky using astrometric information in addition to photometry. GD-1 was immediately detected (Malhan et al. 2018b) as one of the highest contrast stellar streams in the Galactic halo. Additional stars surrounding the stream were detected (Price-Whelan & Bonaca 2018) including an off-track "spur" (marked in Figure 1a), these features may be revealing the effect of massive perturbers (Bonaca et al. 2019) or they may point to the possibility that the progenitor of GD-1 originated within a larger system (Malhan et al. 2019).

THE STREAMFINDER GD-1 SAMPLE

We aim to derive a sample of GD-1 stars from which we will be able to examine its stellar number density profile. Constructing such a sample is not entirely straightforward though, because of the substantial contamination from normal Galactic field populations over the large area of sky that this structure covers. To extract a clean sample, we first need to know the large-scale behavior of GD-1 in position, parallax, proper motion and photometry.

We will therefore begin our analysis by first deriving the properties of GD-1 from a sample of 811 candidate member stars identified with the STREAMFINDER algorithm (Malhan & Ibata 2018;Malhan et al. 2018b). This software provides a means to assign a likelihood to every star in a dataset according to the possibility of whether the star can be grouped with other stars into a stream-like structure. The adopted algorithm parameters are stated in Ibata et al. (2019); in particular, we searched for stream stars down to G 0 = 19.5 mag using a stream template of Gaussian width 0.05 kpc, and of length 20 • . Three different stellar populations models from the PARSEC library (Bressan et al. 2012) were used, with age and metallicity: (12.5 Gyr, -2.0), (12.0 Gyr, -2.3), and (10.0 Gyr, -1.7). For every star we adopted the most likely stream solution obtained from one of these three age-metallicity choices. The resulting spatial distribution of the candidate GD-1 members is shown in Figure 1a, displayed in the φ 1 , φ 2 coordinate system of Koposov et al. (2010), where φ 1 corresponds to position on a great circle that is approximately parallel to the GD-1 stream. (For easier comparison to maps in Equatorial coordinates, the φ 1 axis in all figures is displayed such that φ 1 increases towards the left).

In Figure 2 we show the color-magnitude diagram (CMD) of a sub-sample of the STREAMFINDER detections, using photometry extracted from the second data release (DR2) of the Pan-STARRS survey (Chambers et al. 2016). Since the stream displays a substantial distance gradient, we selected the sub-sample to lie between -50 • < φ 1 < -40 • , where the distance is approximately constant. We have chosen to show Pan-STARRS (instead of Gaia) photometry here because of the much smaller uncertainties at the faint end of the CMD. The two PARSEC stellar population models with age and metallicity (12.5 Gyr, -2.0) and (12.0 Gyr, -2.3) can be seen to give a reasonable representation of GD-1, and are consistent with spectroscopic measurements derived from Segue and LAMOST (Malhan & Ibata 2019). We include the more metal-rich model with (10.0 Gyr, -1.7) to represent an extreme upper limit to the CMD properties of GD-1.

As part of an on-going follow-up survey of the stream stars detected with the STREAMFINDER algorithm, we observed 29 GD-1 candidate stars with the high-resolution ESPaDOnS spectrograph at the CFHT, 5 stars with the EFOSC2 spectrograph at the NTT, as well as 2 stars with the high-resolution UVES spectrograph at the VLT. The ESPaDOnS spectra were extracted and wavelength-calibrated with the Libre-esprit software (Donati et al. 1997), and we used the ESOREX pipeline to perform the same task with the EFOSC2 and UVES spectra. The radial velocities of all 35 stars were measured by cross correlation against the radialvelocity standard HD182572 using the "fxcor" command in IRAF. The average uncertainty of the stars observed with ESPaDOnS and UVES is 0.9 km s -1 . An arti- cle presenting the spectroscopic follow-up survey of the STREAMFINDER detections is currently in preparation, and we defer a detailed exposition of the data to that contribution.

We cross-matched the STREAMFINDER sample against public spectroscopic surveys, finding matches with 2 stars in APOGEE-2 (Majewski et al. 2017), 2 stars in the Gaia Radial Velocity Spectrometer (RVS) sample, 43 stars in LAMOST DR5 (Cui et al. 2012) and 91 stars in SDSS/Segue (Yanny et al. 2009). The final velocity sample (including our CFHT/ESPaDOnS and ESO observations) consists of a total of 156 distinct stars out of the sample of 811. For those stars with multiple measurements, we adopted the measurement that possesses the lowest uncertainty. These velocity measurements are displayed in Figure 1b, along with their uncertainties.

The heliocentric radial velocity of the stream can be seen to change smoothly by almost 600 km s -1 over the 95 • that we detect it over. We performed a simple empirical fit to the velocity data v, rejecting those stars with |v -v fit | > 20 km s -1 + 2δv, where δv is the radial velocity uncertainty. The fitted polynomial

v fit = 90.68 φ 3 1 + 204.5 φ 2 1 -254.2 φ 1 -261.5 (1)
(with velocities in km s -1 and φ 1 in radians) is shown with a blue line, and the 117 stars that are retained in the fit are shown in black, while the 39 rejected stars are colored magenta. We deliberately use empirical fits in the present contribution rather than fitting a stream model so as to avoid mismatch biases from errors in the Galactic potential model.

A further empirical fit S(φ 1 ) is made to the φ 2 trend of the 117 velocity-confirmed members. We find:

S(φ 1 ) = 0.008367 φ 3 1 -0.05332 φ 2 1 -0.07739 φ 1 -0.02007 , (2)
where all angles are in radians. This fit is shown with the solid blue line Figure 1a. The majority of the velocityconfirmed members lie within 0 • .6 of this fitted line (i.e., between the dotted lines): discounting the "spur" feature, only 4 velocity members extend beyond 0 • .6. In contrast, of the 39 velocity non-members, 26 lie beyond 0 • .6. Thus the contamination of the velocity sample within |S(φ 1 )| < 0 • .6 is only 11%, which motivates our choice of selecting stars from the STREAMFINDER sample from within this region of sky.

The proper motion properties of GD-1 are displayed in Figure 3, derived from the 603 STREAMFINDER candidates with |φ 2 -S(φ 1 )| < 0 • .6 and that are not radial velocity outliers. The fitted polynomial relations are:

µ α,fit = 3.794φ 3 1 + 9.467φ 2 1 + 1.615φ 1 -7.844 (3)
and

µ δ,fit = -1.225φ 3 1 + 8.313φ 2 1 + 18.68φ 1 -3.95 , (4)
with φ 1 in radians and the proper motions in mas/yr. The STREAMFINDER software returns the most likely distance solution for the stream model at the position of every star in the sample. The search for stellar streams was initially undertaken using only Gaia photometry, and we conducted our spectroscopic follow-up survey based on those data. However, we have recently updated the software to allow us to include other photometric catalogs. The only conceptual change to the software that this entails is an additional factor in the probability density model of the stream (Equation 2of Ibata et al. 2019) model predictions with a simple Normal distribution, i.e., P color,PS1 = N (x), where x ≡ ((g -r) 0 -(gr) 0,model )/δ(g -r), where δ(g -r) is the color uncertainty. This upgrade to the software significantly decreases the uncertainties on the distance estimates. The resulting distance trend is shown in Figure 4, which we have fit (red line) with the following polynomial:

D(φ 1 ) = -4.302φ 5 1 -11.54φ 4 1 -7.161φ 3 + 5.985φ 2 1 + 8.595φ 1 + 10.36 , (5)
where D is in kpc and φ 1 is in radians. This fit was made to the velocity-confirmed stars (filled black circles), but it clearly also encapsulates the trend of the full STREAMFINDER sample (blue points).

The analysis described so far in this section has allowed us to derive empirical fits to the track of the stream on the sky, to its line of sight velocity profile, to the proper motion gradient in µ α and µ δ , and to the distance gradient. With these ingredients we can now return to the original Gaia catalog and examine the density distribution along the stream.

GD-1 DENSITY PROFILE

While the STREAMFINDER provides a sample of stream members, the reader may be concerned that the algorithm's parameters could bias the results in a compli- cated way. For this reason, we now proceed to extract two additional stream samples selected in a more traditional way, to serve as comparisons.

We first extract a sample of GD-1 stars from the Gaia DR2 catalog, taking those stars with G 0 < 20 mag, that have a full 5-component astrometric solution, and that possess a flux excess E(≡ phot bp rp excess factor) in the range: Lindegren et al. 2018 for the motivation for this constraint).

1 + 0.015(G BP -G RP ) 2 < E < 1.3 + 0.06(G BP -G RP ) 2 (6) (see
For every star in the Gaia DR2 catalog in the region -100 • < φ 1 < 20 • and -10 • < φ 2 < 10 • , we calculate the probability of the star belonging to an idealised stream model, which is simply the empirical sky position, distance, and proper motion profiles convolved with appropriate Gaussians: P stream = P width × P × P color × P color,PS1 × P µ , (7) where P width is the probability that the star is located at the observed φ 2 -S(φ 1 ) perpendicular distance from the stream track, P is the probability of the observed parallax given the distance model, P color is the probability of the observed Gaia G BP -G RP color, given the distance model and the stellar population model, and P color,PS1 is the same probability for the Pan-STARRS photometry. These four probability terms are modelled as one-dimensional Gaussians, and the observed uncertainties are taken into account by adding the uncertainty in quadrature with the intrinsic model dispersion (for the P term, we also adopt the parallax zero-point of -0.029 mas found by Lindegren et al. 2018). The fifth factor P µ is the probability of stream membership given the measured proper motion differences from the model (∆ µα and ∆ µ δ ), and is given by:

P µ = 1 2πσ µα σ µ δ 1 -ρ 2 × exp(- 1 2(1 -ρ 2 ) ∆ 2 µα σ 2 µα + ∆ 2 µ δ σ 2 µ δ - 2ρ∆ µα ∆ µ δ σ µα σ µ δ ) . (8)
We thus take into account the proper motion uncertainties σ µα and σ µ δ their correlation C ≡ pmra pmdec corr (see Lindegren et al. 2018), which is incorporated into the term

ρ = C σ µα σ µ δ (σ 2 µα + w 2 µ )(σ 2 µ δ + w 2 µ ) , (9)
which can be derived by convolving the two-dimensional covariance matrix with an isotropic two-dimensional Gaussian of dispersion w µ .

For the present search, we assume a Gaussian width of the stream of 50 pc, and we (generously) allow a dispersion in proper motion w µ equivalent to 10 km s -1 in velocity. We adopt the same three PARSEC stellar populations models as used above for the STREAMFINDER, with age and metallicity (12.5 Gyr, -2.0), (12.0 Gyr, -2.3), (10.0 Gyr, -1.7), and for each star, we select the solution that yields the highest probability.

The result of the search is shown in Figure 5. In panel (a) we show the 4784 stars that lie within 3σ of the model in parallax, proper motion, and Gaia and Pan-STARRS photometry. While GD-1 is clearly visible, there is a non-negligible amount of contamination in the map, especially towards the extremities of the structure, where it approaches regions of low Galactic latitude. Panel (b) also shows a 3σ membership cut, but this time the distance to the stream track is also taken into account, yielding a sample of 868 stars. A comparison to (a) shows that this 3σ cut corresponds to the sample that one would select visually as probable GD-1 members. Note that this 3σ cut removes the "spur" feature (Price-Whelan & Bonaca 2018) visible in Figure 1. Since we are attempting to ascertain the reliability of the density profile, it is useful to consider the properties of alternative GD-1 samples. Therefore, in panel (c), we present the map of the 1171 stars that lie within 3σ of the empirical model, but this time ignoring the Pan-STARRS photometry. As expected, the contamination is higher in this case.

Figure 6 condenses this information into onedimensional star-count profiles, where the 3σ sample derived with Gaia and Pan-STARRS information is presented in (a), while (b) ignores the Pan-STARRS colors. The estimated contamination in each sample is shown in the blue histograms, and has been subtracted from the density profiles of interest (red histograms). This contamination is estimated by selecting those stars with |φ 2 -S(φ 1 )| > 2 • and with |φ 2 | < 10 • . In (c) we show the profile of the STREAMFINDER sample of 603 stars within |φ 2 -S(φ 1 )| < 0 • .6 (and that are not radial velocity outliers). A simple quadratic was fitted to each profile (dotted lines); these are, respectively, for panels (a,b,c):

C Gaia+PS1 (φ 1 ) = -37.51φ 2 1 -46.51φ 1 + 14.37 C Gaia (φ 1 ) = -50.48φ 2 1 -63.37φ 1 + 13.86 C STREAMFINDER (φ 1 ) = -28.7φ 2 1 -35.96φ 1 + 11.26 , (10)
where φ 1 is in radians and the counts C are per bin of width 2 • .5. Interestingly, the density distribution along the stream displays prominent spikes that can be seen as high contrast peaks above the low-order fit. In Figure 7 we reproduce the sky distribution of these sources in the STREAMFINDER sample, colored according to local density.

It is not surprising that the three distributions in Figure 6a-c are not identical given the different selection procedures. However, they all display at least four extremely prominent peaks at the same locations, and as such they paint a consistent picture of the large-scale properties of the GD-1 system.

DENSITY ARTIFACTS DUE TO MIS-CLASSIFICATION AND INCOMPLETENESS

Before analysing the implications of our measurement of the density profile along GD-1, we will first discuss the limitations of stellar density maps derived from groundbased imaging. Modern wide-field cameras allow one to detect structures of very low surface brightness simply by counting individual resolved stars. Typically, these sources are revealed as small enhancements over the foreground and background contaminating populations. Generally, the image quality of a camera degrades away from the field center, causing both higher photometric uncertainties, and poorer classification constraints, so that the fraction of stars that may be confused with galaxies (and vice-versa) worsens towards the edges of the field of view.

The weather conditions obviously also change over the course of a large survey, leading to varying survey depths as the transparency of the sky changes, as well as different depths for accurate star/galaxy classification. Temperature variations will also lead to variations in the quality of the focus.

All these factors affect the spatial homogeneity of a survey in a complicated way that is not easy to estimate or correct for. This is especially the case for public surveys where the information about the observing conditions that went into producing the data in a particular region of sky are difficult to recover. We therefore felt that it would be useful to investigate how reliably a ground-based survey such as Pan-STARRS could be used to measure large-scale stellar density. We stress We have intentionally omitted the Poisson uncertainties on these distributions, so as to highlight the fact that they are not independent samples. In principle the distributions should be identical; but they are not due to differences in completeness and star/galaxy classification between the two surveys. The substantial (∼ 20%) variations will result in large spatially-dependent errors in any density map derived from such data. This problem is further compounded if the observed counts are dominated by a large contaminating population that needs to be subtracted off (as is the case with GD-1). b: deviations between the counts of point-sources selected within the blue rectangle in Figure 8.

that Pan-STARRS photometry is known to be photometrically extremely well calibrated (with a reliability of 7-12 millimags; Chambers et al. 2016); the issue we wish to assess here is its homogeneity to classification and completeness over large fields.

To this end, we decided to compare the GD-1 stream region in Pan-STARRS to a deeper survey, taken in good seeing conditions with CFHT/MegaCam, which was previously analysed by de Boer et al. (2018). We retrieved the images from the CFHT archive and processed them with the same procedure as applied to data from the Canada-France Imaging Survey (Ibata et al. 2017). The dataset consists of 528 g-band images and 516 r-band images, all of exposure time 50 s, that cover the gray-shaded region in Figure 1a.

The CFHT/MegaCam images were recalibrated onto the Gaia DR2 astrometric reference, which was also used as the astrometric reference for the Pan-STARRS DR2 catalog. The zero-points of the CFHT/MegaCam g and r-band photometry were calibrated onto the Pan-STARRS DR2 survey, adopting the color transformations1 :

g CFHT = g PS1 + 0.014 + 0.059x -0.00313x 2 -0.00178x 3 r CFHT = r PS1 + 0.003 -0.050x + 0.0125x 2 -0.00699x 3 (11)
where x ≡ (g -i) PS1 . The Cambridge Astronomical Survey Unit (CASU) software (Irwin & Lewis 2001) was used to measure the photometry and perform the star/galaxy classification.

Figure 8 shows (in blue) the resulting color-magnitude diagram (CMD) of all stellar sources identified in the CFHT/MegaCam survey. The black points mark the positions of the STREAMFINDER GD-1 candidates that are present in the CFHT/MegaCam survey region; the main sequence turnoff of GD-1 is clearly visible, and the sample also contains some sub-giant and red giant branch stars.

For our comparison test, we decided to isolate the stars in the red selection box in Figure 8, as this corresponds to the CMD location where GD-1 has its highest contrast over the contaminating populations of the Milky Way, and it will be the signal in this CMD region that a matched filter will enhance. A rectangular box is chosen for simplicity, selecting stars with (g -r) 0 ∈ [0.19, 0.24] and g 0 ∈ [18, 19.5]. The average photometric uncertainties of the stars in this box are below 0.01 mag in both surveys and in both colors. A sample of point-sources is selected from the Pan-STARRS DR2 survey by (conservatively) retaining only those sources where the r-band PSF magnitudes agree with the aperture magnitudes to within 0.05 mag. The Pan-STARRS targets are further required to have a minimum of two detections in the g-and r-bands, and to have qualityFlag = 4 (which identifies good-quality measurements in Pan-STARRS).

Figure 9a which contains halo main sequence turn-off stars. Given that the halo is expected to be an ancient, dynamically well-mixed population, this sample should be spatially smooth. Within this box the typical uncertainties of the CFHT/MegaCam and Pan-STARRS photometry are 0.01 mag and 0.05 mag, respectively. Significant differences in the number-counts profiles between the CFHT/MegaCam and Pan-STARRS results are seen again, and a comparison between panels a and b of Figure 9 demonstrates that the deviations do not match up spatially between the samples.

In addition to the large-scale variations of the type seen in Figure 9a, which may be due to variable transparency and seeing over the course of a survey, periodic camera-sized density variations are often seen in widefield maps. Such artifacts can sometimes be spotted following the survey tiling pattern (see, e.g., the ripples in the PAndAS survey in Figure 11 of Ibata et al. 2007, or in the u-band of the SDSS in Figure 3 of Ibata et al. 2017).

The astrophysical interpretation of density variations measured from ground-based wide-field surveys therefore requires a very careful correction for spurious signals.

Of course space-missions may also have spatiallydependent artifacts. In the case of Gaia DR2 there are particularly noticeable stripes of incompleteness that follow the scanning pattern (see e.g. Lindegren et al. 2018), and these artifacts will contribute to the measured spatial variations in density. These problems are difficult to perceive in our previous maps, because of the low density of sources in the stream. However, by examining the spatial distribution of all Gaia sources with G 0 < 20 mag (Figure 10) sufficient statistics are attained to reveal numerous track-like diagonal underdensities crossing the path of GD-1. These are particularly noticeable in the interval φ 1 = [-60 • , -40 • ], where they cause narrow (∼ 0 • .2) dips of ∼ 50% lower density with a periodicity in φ 1 of slightly over 1 • . A wider (∼ 3 • .6) band of lower density is also visible intersecting the stream path at φ 1 = -56 • .5. These artifacts cause a patchy incompleteness in the stream survey, and will contribute spurious gap-like information to the density power spectrum of GD-1 derived from Gaia DR2 data.

MODELLING THE DENSITY PROFILE

Having presented some of our concerns on the limitations of star-counts measurements, we now proceed to model the density profiles measured in Section 4. The star distribution of Figure 7 is strikingly reminiscent of the epicyclic overdensities seen in simulations of slowly disrupting clusters: see especially Figure 7 (2012). This suggests that we can constrain the mass of the progenitor from the periodicity of the density peaks. We consider here the simplest case, in which a satellite moves on a circular orbit, and material is lost at the escape radius with null velocity. Given that GD-1 lies on a low eccentricity orbit (e = 0.33, Willett et al. 2009), this simple configuration is not too unrealistic.

In this situation the distance between two overdensities due to the epicycles along the stream is (Küpper et al. 2008)

y C = - 4πΩ κ 1 - 4Ω 2 κ 2 x E , (12)
where Ω and κ are the circular and epicyclic frequency at the Galactocentric distance r of the cluster and x E is the escape radius from the cluster. The escape radius can be approximated by 2.88 times the Jacobi radius (Varghese et al. 2011;Fardal et al. 2015)

x E = 2.88 × GM 4Ω 2 -κ 2 1/3 , ( 13
)
where M is the total mass of the cluster.

If we assume a spherical Milky Way model with a constant circular velocity v 0 (i.e. with a logarithmic potential), then Ω(r) = v 0 /r, and κ(r) = √ 2Ω(r). Eqs. ( 12 and ( 13) then simplify to

y C = 4π √ 2 x E , (14)
and

x E = 2.88 × GM r 2 2v 2 0 1/3 , (15)
respectively. Assuming that we are observing a distance between the peaks along the stream of y obs , this implies a relationship between v 0 and

M M y 3 obs v 2 0 270.26 × Gπ 3 r 2 . ( 16
)
Inspection of Figure 4 reveals that the distance of the stream in the zone of interest is d ∼ 8 kpc, and the angular distance of the peaks ∆φ ∼ 15 • . This would suggest that y obs = d∆φ ∼ 2 kpc, but this estimate ignores the fact that the stream is not perpendicular to the line of sight. After correcting for the projection effect, we will measure y obs = 2.64 ± 0.18 kpc in Section 7. Given that the stream lies at a Galactocentric distance r ∼ 15 kpc, we then expect the relationship between circular velocity and progenitor mass shown in Figure 11a.

In Figure 11b, we show the separation of the density peaks of GD-1, predicted from Eq. ( 12), taking a circular velocity curve with v 0 = 229 km s -1 (consistent with the measurement of v 0 = 229.0 ± 0.2 km s -1 at the Solar radius by Eilers et al. 2019), and using a distance of GD-1 from the Sun of d = 8 kpc.

Having established plausible masses for the GD-1 progenitor, we now examine whether a disrupting N-body model can give rise to the observed stream density profile. For this N-body simulation, we adopted the Galactic potential of Dehnen & Binney (1998) (their model 1) for the bulge, thin disk, thick disk and interstellar medium. For the dark matter halo, we used a Navarro et al. (1997) model similar to the dark matter halo found recently by Cautun et al. (2019), with a virial radius of 206 kpc, a concentration of c = 12, but with an oblateness of q = 0.82 (as derived by Malhan & Ibata 2019 from modelling GD-1). These choices lead to a dark halo mass of 9.6 × 10 11 M . With this Galactic potential model, the circular velocity at the Solar radius (R = 8.129 kpc, Gravity Collaboration et al. 2018) is 229 km s -1 . We integrated backwards in time for 2 Gyr starting from (R.A., Dec) = (157 • .6, 43 • .71667), d = 8.25 kpc, (µ α , µ δ ) = (-6.53 mas/yr, -11.0 mas/yr), and v helio = -90 km s -1 . We then integrated a Kingmodel (King 1966) forwards in time for 2 Gyr, using the gyrfalcON N-body integrator (Dehnen 2000).

The King model was set up to produce a rapidlydisrupting structure so that at the end-point of the simulation there would be no discernible bound structure. The initial mass of the best model we found is 3 × 10 4 M (a factor of ∼ 3 lower than our prediction from Figure 11). The model also possesses a central potential W = 3.0 and a King model tidal radius of r t = 0.17 kpc. We used 50, 000 particles, and a softening length of 1.5 pc.

In order to account for the incompleteness of the Gaia DR2 survey, we applied the Gaia completeness "flat field" map (shown in Figure 10b, additionally excising a 2 • circle around M67) to the final simulation output.

At the end of the simulation the spatial structure of the stream follows the derived large-scale threedimensional properties of GD-1 fairly well, as we show in Figure 12. Although the cluster is completely disrupted by this point (as observed), we estimate that the position of the progenitor in this model if it had survived would have been φ 1 = -29 • .925, φ 2 = 0 • .096. Several strong peaks can be seen to be present in the density profile, as shown in Figure 6d. Such peaks are a generic property of simulated streams that dissolve slowly in this way (Küpper et al. 2012).

Our limited exploration of the parameter space of the simulations suggests to us that it is challenging to match an N-body stream model to these observations, in part because of the cubic dependency of the tidal radius on the progenitor's mass (Eq. 13), and because of the rapid time evolution from what must have been a bound structure to complete dissolution. This renders the location and contrast of the peaks in the N-body simulation very sensitive to the modelled initial conditions.

Despite displaying a multi-peaked density profile, the best N-body model we have found so far does not reproduce faithfully the observed peak morphology (Figure 6d). It is particularly noticeable that the peaks are wider than in reality, and the narrow peak at φ 1 ∼ 0 • is not present. Some of our N-body models do produce the density spike at φ 1 ∼ 0 • , but obtaining that peak comes at the cost of much lower peak contrast elsewhere in the profile. We are currently in the process of simulating a large library of such models, which will be presented in a future contribution. It is likely that the slight discrepancies between the path of the stream through the Galaxy in the simulation and in the observations that can be seen in Figure 12 are due to the adopted Galactic potential model giving a slightly incorrect acceleration field; this also will be explored in future work.

We suspect that the large (factor of ∼ 3) overestimate of the progenitor's mass made by the analytic model (Eq. 15) compared to our best N-body simulation is because the assumptions underlying that model do not hold true. In particular, the assumption of constant mass is obviously a poor one in relation to a structure that ends up disintegrating completely.

POWER SPECTRUM ANALYSIS

For completeness, we finally calculate the power spectrum of the tidal stream following Banik et al. (2019). However, we feel that at present the sky position of the remnant of the GD-1 progenitor remains highly conjectural (Price-Whelan & Bonaca 2018; Malhan et al. 2018a;de Boer et al. 2018de Boer et al. , 2019;;Webb & Bovy 2019), so splitting the structure into leading and trailing arms is not justified. We use the csd algorithm in scipy to calculate the density power spectrum of the profiles shown in Figure 6a-c normalised by the respective quadratic fits to the continuum. The result of this calculation is shown in Figure 13, as a function of inverse wavenumber 1/k φ1 in the φ 1 coordinate. The uncertainties on the power spectrum are derived by re-running the procedure on 1000 randomly-drawn profiles consistent with the profile uncertainties. The blue and green lines show the power spectra of the Gaia+Pan-STARRS and Gaiaonly profiles (from Figure 6a and 6b, respectively), while the yellow line is the power spectrum derived from the STREAMFINDER sample. The three samples show similar, but not identical, behavior.

To serve as a comparison, we also calculate the power spectrum of the Galactic halo contamination in this region of sky (red line). For this, we chose to use the contamination profile previously shown in Figure 6a (blue histogram), derived from the sky region with |φ 2 -S(φ 1 )| > 2 • which contains 861 stars in the range φ 1 = [-90, 10], an almost identical number to the Gaia+Pan-STARRS GD-1 sample (868 stars). The power spectrum of the smooth halo contaminants can be seen to be very similar to that of the STREAMFINDER sample.

The measured power at different angular scales may come from clear structures, such as the spikes seen in Figure 6, and some signal may be due to interactions with invisible dark matter sub-halos. However, it is also possible that Gaia's scanning law (Figure 10) has imprinted a signal on the power spectra, and the fact that the background shows a similar density power spectrum to GD-1 lends weight to this concern. Indeed, the power spectra of the background and STREAMFINDER samples fully overlap within 1σ over all spatial scales probed. Note in particular that the background sample extends over a very much wider range in the φ 2 coordinate than the GD-1 stream samples, which was necessary in order to extract a similar number of stars. This means that the imprint of the scanning law will probably be diminished in the background sample since any artifacts will be averaged-out over the large φ 2 range.

In Figure 13 we also show the power spectrum of the N-body model presented in Section 6. Despite the fact that the model has been integrated in a perfectly smooth Galactic potential, the power spectrum displays a considerable similarity to the observed profiles.

The slight bump in the power spectra at ∼ 15 • in Figure 13a is due to the strong peaks seen in Figure 6. Because of projection effects, the periodicity of the fea- . Density power spectra of the stream density profiles from Figure 6. These are calculated using the Welch (1967) method, as provided by the csd algorithm in scipy. The (1σ) uncertainties on each power spectrum have been estimated by re-sampling the corresponding density profiles 1000 times. For clarity, the power spectra have been separated into two panels, and the results for the STREAMFINDER sample are reproduced in both panels to allow easier comparison. The power spectrum of the background sample (red) corresponds to the contamination profile estimated for the Gaia+PS1 sample. The black line is derived for the N-body model presented in Section 6.

tures is somewhat veiled when they are examined in the φ 1 coordinate, but it becomes obvious after changing coordinates to a proper path length along the stream. Defining

s(φ 1 ) = φ1 0 D 2 + dD dφ 1 2 dφ 1 (17)
to measure this path length, we re-calculate the density power spectra as a function of inverse wavenumber 1/k s in this s coordinate, and show the results in Figure 14. The epicyclic over-densities produce a very clear signal at 1/k s = 2.64 ± 0.18 kpc (calculated by fitting a Gaussian to the 1000 random realisations of the STREAMFINDER sample, selecting data in the range 1/k s = [1.3, 4.5] kpc), and all three samples shown in Figure 14a are consistent with each other. In Figure 14b, we now see that the background sample does differ from the STREAMFINDER sample, but mostly by the fact that it does not exhibit a pronounced peak at 1/k s = 2.64 ± 0.18 kpc. At small separations (where Banik et al. 2019 find that the perturbing influence of LCDM substructures is required), the background and STREAMFINDER samples have identical behavior.

To interpret these power spectra it is useful to know the level of the shot noise. To estimate this, we made 1000 realisations of a uniform distribution in s containing 868 stars, and calculated the corresponding power spectra. The resulting distributions are flat in 1/k s and possess mean and 1σ uncertainties as displayed in Figure 14a (black errorbar).

We note in passing that the larger epicyclic peak distance in our N-body model (black line in Figure 14b) suggests that the N-body model was ∼ 30% too massive at the time when the peaks were formed (estimated from Eq. 16). This hint will be explored in future work.

CONCLUSIONS

We have used the STREAMFINDER algorithm to isolate a sample of stars of the GD-1 stellar stream over ∼ 100 • of the Northern sky. Our radial-velocity follow-up of these STREAMFINDER candidate members shows that there is very little contamination (∼ 10%) if the sample is spatially restricted to being close (< 0.6 • ) to the fitted path of GD-1. We take advantage of this sample to fit empirical relations to the sky position, radial velocity, proper motion and distance to the structure, which are then used to extract two other clearly-defined samples of GD-1 stars from the Gaia DR2 and Pan-STARRS catalogs. As Figure 13, but as a function of distance s along the stream, rather than observed angle φ1. The fact that the peak at 2.64 ± 0.18 kpc has much higher contrast here implies that projection effects smear out the periodic signal in Figure 13. The black errorbar in (a) shows the 1σ uncertainty level due to shot noise in an artificial stream with 868 stars uniformly distributed in s.

We find that the three different GD-1 samples we have constructed have a similar spatial distribution. In particular, very strong peaks are present along the stream, spaced by 2.64 ± 0.18 kpc, with a contrast exceeding 3:1.

While the density power spectrum may in principle contain information about the prevalence of perturbing massive bodies, for spatial separations up to 2 kpc, we find a very similar behavior between the GD-1 samples we examined and a background profile extracted from the sky regions immediately adjacent to GD-1.

We also present a comparison between star-counts profiles derived from Pan-STARRS DR2 and a deeper survey with the MegaCam wide-field camera at the Canada-France Hawaii Telescope taken in much better seeing conditions. In a color-magnitude region where the GD-1 stream has highest contrast over the contaminants, we find substantial spatially-dependent differences in the corresponding star counts. We attribute these differences to variations in the observing conditions, leading to variations in completeness and variations in star/galaxy discrimination between the two imaging surveys. These errors are very hard to identify and correct for in large ground-based surveys without an external deeper dataset, and may strongly affect conclusions of the prevalence of density gaps in streams.

In contrast, space-based surveys may be more powerful for measuring density profiles, because they are unaffected by our variable weather. However, with a scanning instrument like Gaia, there may be spatially varying incompleteness due to the way in which the survey has been designed to cover the sky. In the particular region around GD-1 (Figure 10) there is evidently significant incompleteness on a range of spatial scales.

Nevertheless, the strong peaks detected here (Figures 6 and 7) are clearly real, being also visible in the matched filter maps presented in the discovery paper (Grillmair & Dionatos 2006). Our modelling shows that these features are most probably due to epicyclic motion in the stream. This conclusion is strengthened by the finding of a strong periodic signal in the density power spectra when we correct for the projection effects (Figure 14). The fact that these periodic density variations are still visible, and have not yet been washed out (as stars mix over time due to the dynamical evolution of the stream) implies that the progenitor of the system went through its final disruption stage only very recently.

In order to obtain reliable constraints of the effect of dark matter substructures from the density profiles of GD-1 it will be necessary to fully account for these internal dynamical properties of the stream, as well as the external perturbations from the baryonic components of the Milky Way. Finally, the (probably very complicated) instrumental sensitivity function will need to be corrected for.

Figure 1 Figure 2 .

 12 Figure1. a: Spatial distribution of GD-1 stream star detections, using the STREAMFINDER software. The φ1, φ2 coordinate system ofKoposov et al. (2010) is used, where φ1 points approximately along the stream, while φ2 is perpendicular to it. The points are color-coded according to the AV extinction, which can be seen to change in a complex manner along this long structure, possessing wave-like variations on scales of degrees. Additionally, we mark (in gray) the CFHT/MegaCam survey region discussed in Section 5. The large foreground open cluster Messier 67 strongly contaminates the GD-1 stream in the marked circular region (which appears elliptical due to the stretching of the φ2 axis). The GD-1 "spur" structure (Price-Whelan & Bonaca 2018) is also marked. b: Stars with radial velocities measurements. The magenta open circles show the GD-1 candidates identified by the STREAMFINDER, whose measured velocities are incompatible with being GD-1 members. The probable GD-1 stars are displayed in black. These larger open circles (with the same color-coding) are also shown in panel (a). Note that the velocity outliers have a strong tendency to also be spatial outliers. Furthermore, the probable members (based on their radial velocities) define a narrow spatial sequence, with the exception of the "spur" grouping. (The error bars show 1σ velocity uncertainties).

Figure 3 .

 3 Figure 3. Proper motion profiles in µα (a) and µ δ (b). The STREAMFINDER sample shown here (containing 603 stars) has been trimmed spatially to lie between the dotted lines in Figure 1a, where contamination is low. The fitted cubic polynomials are defined in the text. (1σ uncertainties are displayed).

  photometric distance D [kpc]

Figure 4 .

 4 Figure 4. Heliocentric distance profile along GD-1, derived by the STREAMFINDER using Pan-STARRS photometry in addition to Gaia DR2 data. The red line shows a quintic polynomial fit to the velocity-confirmed data (black dots). (Error bars show 1σ distance uncertainties, as estimated by the software).

Figure 5 .

 5 Figure 5. Selection of GD-1 candidates directly from the Gaia DR2 and Pan-STARRS DR2 catalogs according to their proximity to the empirical stream model. (a) shows the sources within 3σ of the model, but ignoring the sky proximity criterion (P width ). Moderate contamination can be seen towards the ends of the stream. (b) as (a) but now accounting for the proximity to the S(φ1) track. Using only Gaia photometry (c) leads to somewhat higher contamination.

Figure 6 .

 6 Figure 6. Counts (with 1σ uncertainties) of candidate members as a function of φ1. (a) shows the density profile based on Gaia data, complemented with Pan-STARRS photometry, based on the map in Figure 5a. The blue histogram shows the expected contamination level derived from the sky region at |φ2 -S(φ1)| > 2 • . The quadratic fit (dotted line) is defined in the text. (b) shows the same information as (a), but ignoring the PS1 information. (c) displays the profile derived from the STREAMFINDER detections with |φ2 -S(φ1)| < 0 • .6, where the radial velocity non-member stars identified in Figure 1 have been rejected. (d) shows the profile of the N-body simulation presented in Section 6, which possesses some of the main features seen in the observations.

Figure 7 .

 7 Figure 7. Map of the STREAMFINDER GD-1 sources color-coded by local density enhancement over a quadratic fit to the counts shown in Figure 6c.

Figure 8 .

 8 Figure8. Color-magnitude distribution of sources in the CFHT/MegaCam survey (blue dots). The stars identified as candidate members in this region by the STREAMFINDER algorithm are highlighted with larger black dots. These can be seen to follow a color-magnitude sequence that includes the main sequence turnoff (red rectangle), where the GD-1 population has the highest contrast over the Galactic contamination. The blue rectangle selects a small dense portion of the main sequence turn-off of the halo.

Figure 9 .

 9 Figure 9. a: Counts as a function of φ1 of stars within the red selection box of Figure 8. The CFHT/MegaCam sample (green) is compared to the Pan-STARRS sample (blue), both selected with identical color-magnitude criteria (the Pan-STARRS bands are converted to CFHT colors, as explained in the text).We have intentionally omitted the Poisson uncertainties on these distributions, so as to highlight the fact that they are not independent samples. In principle the distributions should be identical; but they are not due to differences in completeness and star/galaxy classification between the two surveys. The substantial (∼ 20%) variations will result in large spatially-dependent errors in any density map derived from such data. This problem is further compounded if the observed counts are dominated by a large contaminating population that needs to be subtracted off (as is the case with GD-1). b: deviations between the counts of point-sources selected within the blue rectangle in Figure8.

  shows a comparison of the counts in the CMD selection box along the length of the stream section where the CFHT/MegaCam imaging was obtained. Substantial ∼ 20% variations are seen in what one might have supposed to have been almost identical overlapping samples. Figure 9b repeats this test, but for a slightly fainter sample (selected within the blue rectangle in Figure 8 with (g -r) 0 ∈ [0.25, 0.35] and g 0 ∈ [20, 21]),

  of Küpper

Figure 10 .Figure 11 .

 1011 Figure 10. Spatial artifacts in Gaia DR2 in the vicinity of GD-1. Panel (a) shows the density of all stars to G0 = 20 mag in pixels of 0 • .1 × 0 • .1. The corresponding count level is displayed below the image. The diagonal dark stripes are regions with a lower number of observations due to Gaia's scanning pattern. In (b) we present a "flat-field" map, constructed by dividing the distribution in (a) by a smoothed version of itself (using a Gaussian kernel with σ = 1 • ). The track of GD-1 in these coordinates (blue) is reproduced from Figure 1a.

Figure 12 .

 12 Figure 12. Positions of N-body model particles (orange) are compared to positions of the GD-1 stars (blue), as as derived using the STREAMFINDER distances shown in Figure 4.

  Figure13. Density power spectra of the stream density profiles from Figure6. These are calculated using theWelch (1967) method, as provided by the csd algorithm in scipy. The (1σ) uncertainties on each power spectrum have been estimated by re-sampling the corresponding density profiles 1000 times. For clarity, the power spectra have been separated into two panels, and the results for the STREAMFINDER sample are reproduced in both panels to allow easier comparison. The power spectrum of the background sample (red) corresponds to the contamination profile estimated for the Gaia+PS1 sample. The black line is derived for the N-body model presented in Section 6.

  Figure14. As Figure13, but as a function of distance s along the stream, rather than observed angle φ1. The fact that the peak at 2.64 ± 0.18 kpc has much higher contrast here implies that projection effects smear out the periodic signal in Figure13. The black errorbar in (a) shows the 1σ uncertainty level due to shot noise in an artificial stream with 868 stars uniformly distributed in s.

see http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/megapipe/ docs/filt.html

ACKNOWLEDGMENTS

We would like to thank the anonymous referee for several very insightful comments and suggestions. RI, BF, NM and GM acknowledge funding from the Agence Nationale de la Recherche (ANR project ANR-18-CE31-0006, ANR-18-CE31-0017 and ANR-19-CE31-0017), from CNRS/INSU through the Programme National Galaxies et Cosmologie, and from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No. 834148). Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council of Canada, the Institut National des Sciences de lUnivers of the Centre National de la Recherche Scientique of France, and the University of Hawaii. Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Science de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. Based on observations collected at the European Southern Observatory under ESO programme 103.B-0568(A) and 103.B-0568(B). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. 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/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.