CONSTRAINING THE MILKY WAY POTENTIAL WITH A SIX-DIMENSIONAL PHASE-SPACE MAP OF THE GD-1 STELLAR STREAM ^*

The Sloan Digital Sky Survey (SDSS) is an imaging and spectroscopy survey which mapped a quarter of the sky near the North Galactic Cap. The data have proven extremely useful for the understanding of the Milky Way (MW) halo. In addition to a large list of MW satellites (Belokurov et al. 2007c; Koposov et al. 2007b; Irwin et al. 2007; Walsh et al. 2007), several extended stellar substructures in the MW halo have been found in the SDSS data, such as the tidal tail of the Palomar 5 globular cluster (Odenkirchen et al. 2001; Grillmair & Dionatos 2006a), the Monoceros ring (Newberg et al. 2002), two northern tidal arms of the disrupting Sagittarius galaxy (Belokurov et al. 2006), the so-called Orphan stream (Grillmair 2006; Belokurov et al. 2007b), the Aquila overdensity (Belokurov et al. 2007a), and the very long thin stellar stream called GD-1 (Grillmair & Dionatos 2006b). Recently, Grillmair (2009) claimed the discovery of another four stellar streams. Streams are presumed to be remnants of tidally disrupted satellite galaxies and clusters. They provide important insights into the history of accretion events and the physics of Galaxy formation. The tidal debris from disrupted satellites (clusters) spreads out in an orbital phase on a path that is close to the orbit of the progenitor. Streams tracing out orbits therefore provide opportunities to constrain the MW's gravitational potential.

After initial searches for tidal tails of globular clusters (e.g., Grillmair et al. 1995), it was the extended Sagittarius tidal tail that first made derivation of such constraints practical (see e.g., Ibata et al. 2001; Helmi 2004; Johnston et al. 2005; Law et al. 2005). However, the tidal tail of the Sagittarius galaxy is quite wide and contains a considerable mixture of different stellar orbits, making it complex to model. For constraining the gravitational potential, a stellar stream that is very thin but of large angular extent is ideal because it permits precise orbital models.

The first studies of globular cluster tidal debris only revealed short (≲ 1°) signs of tails, but in recent years with the advent of large photometric surveys such as SDSS and 2MASS and significant advances in the techniques used to find streams, significant progress has been made. The matched filter technique (Odenkirchen et al. 2001; Rockosi et al. 2002) has revealed the beautiful tidal stream of Palomar 5. Detailed analysis of the Pal 5 stream, including kinematics (Odenkirchen et al. 2001, 2003, 2009), has shown the promise of this approach, but also revealed that data over more than 10° on the sky are needed to place good constraints on the potential. Grillmair (2006), Grillmair & Dionatos (2006a, 2006b), and Grillmair & Johnson (2006) were successful in the detection of very long stellar streams using this technique, including the 63° long stellar stream GD-1. Besides the stream length and the approximate distance, most of the properties of GD-1 were unknown. Since the stream is long but relatively thin, with no apparent progenitor remnant, it was suggested that it arose from a globular cluster. In this paper, we make an attempt to determine all possible properties of the GD-1 stream including distance, position on the sky, proper motion, and radial velocity and try to constrain the Milky Way potential using that information. This work goes in parallel with the work done by Willett et al. (2009), but we are able to get a full 6D phase-space map of the stream and to use that map to provide significant constraints on the MW potential. See also Eyre & Binney (2009) for theoretical discussion of using thin streams in order to constrain the MW potential.

In performing this study, we have obtained the first 6D phase-space map for a kinematically cold stellar stream in the MW. We view our present analysis in same sense as a pilot study for the Global Astrometric lnterferometer for Astrophysics (GAIA; Perryman et al. 2001) age, when this ESA space mission will deliver dramatically better data on streams such as GD-1.

This paper is organized as follows. In Section 2, we discuss the analysis of the SDSS photometry, which entails mapping the GD-1 stream in three dimensions (3D) as well as determining its stellar population properties. In Section 3, we present the kinematics with proper motions from SDSS-USNOB1.0 and line-of-sight velocities from SDSS and Calar Alto. In Section 4, we combine this information in an iterative step that involves improved stream membership probabilities, which in turn affects the estimates of proper motions and distances. This procedure results in the most comprehensive six-dimensional (6D) data set for a stellar stream in our MW. In Section 5, we model the stream data by a simple orbit in a simple parameterized gravitational potential. We measure the potential circular velocity and find that the overall MW potential at the GD-1 stream position is somewhat flattened, but that much of that flattening can be attributed to the disk.

The probability that a star is a member of the GD-1 stream depends on its 6D position and its metallicity. In the space of photometric observables, this means that it depends on (α, δ), magnitude, and color. In practice, the determination of the stream's angular position, distance, and metallicity (presuming it is "old") is an iterative process which we detail here.

Grillmair & Dionatos (2006b) made the initial map of the stream using a matched color–magnitude filter based on the color–magnitude diagram (CMD) of M13 observed in the same filters. Not presuming a particular metallicity (e.g., that of M13), we start our analysis with a simple color–magnitude box selection for stars (0.15 < g − r < 0.41 and 18.1 < r < 19.85). The resulting distribution is shown in Figure 1. That particular color–magnitude box was selected as appropriate to find metal-poor main sequence (MS) stars at a distance of ∼ 10 kpc, and indeed the stream is marginally discernible in the figure. With just a color–magnitude box, however, the detection fidelity of that stream is noticeably lower than that achieved by (Grillmair & Dionatos 2006b, their Figure 1). The distribution of stars in Figure 1 is plotted in a rotated-spherical coordinate system (ϕ₁, ϕ₂), approximately aligned with the stream, where ϕ₁ is longitude and ϕ₂ is the latitude. The north pole of that coordinate system is located at α_p = 345987, δ_p = 297331, the zero-point for ϕ₁ is located at α = 200°, and we will use this coordinate system for convenience throughout the paper to describe stream positions (the transformation matrix from (α, δ) to (ϕ₁, ϕ₂) is given in the Appendix).

Zoom In Zoom Out Reset image size

If we integrate the low-contrast two-dimensional (2D) map in Figure 1 along the ϕ₁ axis, creating a one-dimensional profile of the stream, the presence of the stream becomes very clear. Figure 2 shows this profile for stars with 0.15 < g − r < 0.41, 18.1 < r < 19.85, and −60° < ϕ₁ < −10°. In that figure, we also overplot the Gaussian fit to this profile with ∼ 600 stars and Gaussian width ($\sigma _{\phi _2}$) of ∼ 12'. This number of stars corresponds to a total stellar mass of M_* ≈ 2 × 10⁴ M_☉, if we assume a distance of ∼10 kpc (see below), and a Chabrier initial mass function (IMF; Chabrier 2001) with an old, metal-poor stellar population. Given that number of stars, we expect to see around 3000 stream stars in SDSS with r < 22. The mean surface brightness of the stream is around 29 mag arcsec⁻².

Zoom In Zoom Out Reset image size

We expand this approach to the determination of the CMD of the stream. The data and the fit shown in Figure 2 were obtained for a fairly wide color–magnitude selection box. But, we can construct such a profile for any other,e.g., color–magnitude box and that profile can then be fitted by

$$\begin{eqnarray} N_{\rm obs}(\phi _2| {\rm CMD})&=& N_{BG}({\rm CMD})+N_{\rm stream}({\rm CMD})\nonumber \\ &\times & \frac{1}{\sqrt{2\pi }\sigma _{\phi 2}} {\rm exp}\left(-0.5\left(\frac{\phi _2-\phi _{2,0}}{\sigma _{\phi 2}}\right)^2\right),\qquad \end{eqnarray} \tag{ 1 }$$

where CMD refers to a given color–magnitude bin, and where we assume that both center (ϕ_2,0) and width (σ_ϕ2) of the stream are fixed at 0 and 12 arcmin. A fit of the Equation (1) model to the observed data N_obs(ϕ₂|CMD) can be performed in χ² sense. As a result N_stream(CMD), the number of stream stars (and its error), can be determined for each given color–magnitude bin, resulting in a Hess diagram for different pairs of SDSS filters (u − g, g − r, r − i, i − z). Figure 3 shows the resulting Hess diagram of the stream derived in several bands. These clearly show a main sequence (MS). The location of the MS turn-off cannot be clearly identified, although there may be a hint at g = 18.5 and u − g = 1. In Figure 3, we also overplot the Marigo et al. (2008) isochrones with age = 9 Gyr, log(Z/Z_☉) = −1.4 at 8.5 kpc, which seem to match quite well. Ivezić et al. (2008) recently showed that the location of MS stars in the (u − g) − (g − r) color–color plane is a good metallicity diagnostic. Therefore, we construct the (u − g) − (g − r) color–color diagram of the stream stellar population shown in Figure 4, which exhibits a distinct concentration of stars at (0.8, 0.3). This argues for a population of single or a dominant metallicity and we can convert this color location to a metallicity using Equation (4) from Ivezić et al. (2008): [Fe/H]_phot = −1.9 ± 0.1. This provides a metallicity estimate that is directly linked to SDSS spectral metallicity estimate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To derive the metallicity, age, and distance of stream stars in a systematic way, we fit the color–magnitude diagrams using a grid of isochrones populated realistically according to the IMF (Dolphin 2002; de Jong et al. 2008). We focus on fits to the color–magnitude diagrams in u, g and g, r filters, since that the u − g color of the MS turn-off is a good metallicity indicator (Ivezić et al. 2008). We create the synthetic Hess diagrams for a grid of model stellar populations (Girardi et al. 2000; Marigo et al. 2008)⁵ with different ages (3–12 Gyr), metallicities (Z = 0.0001–0.025), distances (6–14 kpc), and a Chabrier IMF(Chabrier 2001). We then explore that grid by computing log-likelihood of the distribution of stars in color–magnitude space. Figure 5 shows the 2D profile likelihoods contours of the age versus metallicity, age versus distance, and distance versus metallicity planes. The filled circle indicates the best-fit model: age = 9 Gyr, log(Z/Z_☉) = −1.4, and distance = 8 kpc. Clearly the age is the least well-constrained parameter; the distance seems to be relatively well-constrained, but has a covariance with [Fe/H]. We will revisit this issue later, as the analysis of Figure 4 implies a lower metallicity. Figure 3 shows that the isochrones are reproducing the observed Hess diagrams well, and hence further in the paper, we will use t = 9 Gyr, log(Z/Z_☉) = −1.4 as the baseline model for the stream's stellar population. It should be noted that the distance measurement from Figure 5 represents the average distance along the stream −50° < ϕ₁ < −20°. In Section 4, we will present estimates of the distance to different parts of the stream.

Zoom In Zoom Out Reset image size

It is noticeable that the metallicity derived from the CMD fitting is higher than from the estimate based on empirical calibration of (Ivezić et al. 2008; see above) and higher than the measurement based on the SDSS spectra given by Willett et al. (2009). This discrepancy is understandable given the known inaccuracies of the isochrones in the SDSS photometric system (An et al. 2008). In particular, Figure 19 of An et al. (2008) clearly shows a mismatch between the fiducial isochrone derived for the M92 globular cluster (which is used elsewhere as a good approximation of old metal-poor stars in the halo) and the theoretical isochrones. For main-sequence stars below the turn-off (which are of the most interest here), the mismatch between a fiducial isochrone of the M53 globular cluster (which has metallicity [Fe/H] ∼−2) from An et al. (2008) and the isochrone which we are using can be approximated by a distance shift of ∼10%. Therefore, in our analysis we reduce all the distances derived on the basis of the CMD fit by 10%. Some remaining systematic error in distances may still exist, due to the described inaccuracies with the isochrones in SDSS filters, although the study of Eyre (2009) seems to indicate it is small. Until isochrones in SDSS filters are fixed, the usage of fiducial isochrones may give more consistent results, but for the current paper we decided to proceed with the theoretical isochrones and 10% distance correction.

Splitting the CMD data into ϕ₁ bins shows that there is a distance gradient along the stream: Figure 6 shows two Hess diagrams obtained for two different pieces of the stream, on the left −40° < ϕ₁ < −20°, and on the right −10° < ϕ₁ < 10°. The baseline isochrone is shifted to distances of 8.5 kpc (left) and 11 kpc (right), respectively. Clearly, the part of the stream at −10° < ϕ₁ < 10° is further from the Sun (as already noted by Grillmair & Dionatos 2006b; Willett et al. 2009).

Zoom In Zoom Out Reset image size

The determination of the CMD properties for the stream allows us to select the possible stream member stars with much less background contamination, compared to a simple color–magnitude box. Figure 7 shows the map of the stream (in ϕ₁, ϕ₂ coordinates) after applying a matched filter based on the CMD of the stream. In Figure 7, each star in the SDSS data set was weighted by the ratio of the stream membership probability and the background probability $\frac{P_{\rm stream}(u-g,g-r,r-i,r)}{P_{BG}(u-g,g-r,r-i,r)}$,where P_stream(u − g, g − r, r − i, r) have been computed based on the stellar population fit (Figure 5), and P_BG(u − g, g − r, r − i, r) have been computed empirically from the regions adjacent to the stream (see, e.g., Rockosi et al. 2002, for the application of similar method). The resulting image after the CMD weighting shows the stream with obviously greater contrast than Figure 1. It is noticeable that the density of stars in the stream does not seem do be constant and even does not change monotonically along the stream, instead the stream seems to have a few clumps and holes. The nature of the substructure in the stream is unclear. It may be related to either the history of the disruption process (Küpper et al. 2010) or, e.g., interaction with dark matter subhalos around MW (Carlberg 2009).

Zoom In Zoom Out Reset image size

In this section, we describe how we derived estimates of the 3D kinematics of the stream by looking at the proper motions and radial velocities of the probable stream members.

3.1. Proper Motion

Despite the distance of ∼ 10 kpc to the GD-1 we demonstrate in this section that it is possible to derive good constraints on its proper motion. The analysis is based on proper motions derived from combining USNO-B1.0 (Monet et al. 2003) with SDSS data (Munn et al. 2004, 2008), which we take from SDSS DR7(Abazajian et al. 2009).

At a distance of ∼ 10 kpc, a stream that is moving perpendicular to the line sight with the velocity of 200 km s⁻¹ (roughly a characteristic velocity in the halo) should have a proper motion of the order of ∼4 mas yr⁻¹. Hence, the expected proper motion is comparable to the precision of individual proper motion measurements of 3–4 mas yr⁻¹ (Munn et al. 2004). As stream member stars that are adjacent on the sky should have (nearly) identical proper motions, we can determine statistical proper motions for ensembles of stars.

We start by deriving the $\vec{\mu }$-distribution of likely member stars, by extending the analysis of Section 2 and using both the angular position on the sky (specifically ϕ₂) and the location in CMDs for each star to evaluate its stream membership probability. Specifically, we grid the proper motion space into "pixels," then we select the stars within each proper motion "pixel" and with high membership probabilities in the CMD (u − g, g − r, r − i, i) space P_stream/P_BG ≳ 0.1. For that sample of stars within each proper motion "pixel" $(\vec{\mu }, \vec{\delta \mu })$, we determine $N_{\rm stream}(\vec{\mu }|\phi _2, {\rm CMD})$ (number of stream stars) by integrating along the stream direction ϕ₁ and fitting the resulting ϕ₂ distribution with the Equation (1). The resulting distribution $N_{\rm obs}(\vec{\mu })$ is shown in Figure 8. The grayscale in the left panel of the figure shows the proper motion distribution of μ_α, μ_δ of probable stream member stars (N_stream) integrated along the stream, while the contours show the proper motion distribution of the background stars selected with the same color–magnitude criteria (corresponding to N_BG from Equation (1)). It is clear that the stream stars are on average moving differently than the bulk of background stars. However, the observed proper motions contain the reflex motion of Sun's motion in the Galaxy. We can account for this and then convert $\vec{\mu }$ to $(\mu _{\phi _1},\mu _{\phi _2})$—the proper motion along the stream $\mu _{\phi _1}\equiv \mu _{\rm along}$ and proper motion across the stream $\mu _{\phi _2}$ in the Galactic rest frame (where ϕ₁, ϕ₂ are the stream coordinates introduced in Section 2). The proper motion component arising from the Sun's movement in the Galaxy can be easily computed for each star.

$$$ \vec{\mu }_{\rm reflex}=\frac{1}{4.74\, \vert \vec{r}\vert } \left(\vec{V}_\odot -(\vec{V}_\odot \cdot \vec{r})\frac{\vec{r}}{ \vert \vec{ r }\vert ^2 }\right), $$$

where $\vec{V}$ is a 3D velocity of the Sun (∼ 220 km s⁻¹) and $\vec{r}$ is the vector from the Sun towards each star. As we approximately know the distance to the stream, this correction $\mu _{\phi _{1,2},c}=\mu _{\alpha,\delta }-\mu _{\rm reflex}$ can be done. We will discuss the consequences of the uncertainties in the Sun's motion and the stream differences in Section 4.

Zoom In Zoom Out Reset image size

The right panel of Figure 8 shows the distribution of $\mu _{\phi _1},\mu _{\phi _2}$ of the stream stars. The contours show the corresponding distribution of background stars with the similar color magnitudes to the stream stars. Reassuringly, we see that stream stars are moving approximately along ϕ₁, i.e., along the stream orbit, an important plausibility check for the correctness of the proper motion measurement. In contrast, the proper motion distribution of the background stars after subtracting the proper motion due to Sun's movement is centered around $(\mu _{\phi _1},\mu _{\phi _2})=(0,0)$, which appears reasonable since with our color–magnitude selection we are selecting primarily the halo stars at distances ∼10 kpc. Those show little net rotation (Carollo et al. 2007; Xue et al. 2008). The estimate $\langle \mu _{\phi _1}\rangle \approx -8$ mas yr⁻¹ also immediately implies that the stream is moving retrograde with respect to the MW's disk rotation.

In Figure 9, we illustrate the proper motion variation along the stream. These plots showing only $\mu _{\phi _1}(\phi _1)$ were derived in the same way as Figure 8, except that we did not integrate in ϕ₁ along the entire stream but only in ϕ₁ intervals. The right panel of Figure 9 shows the distribution of the proper motions along the stream of the background stars. The left panel of the figure shows the distribution of proper motions of likely stream member stars as a function of angle along the stream (the proper motion due to the Sun's motion was subtracted). The left panel reveals a slight but significant gradient in $\langle \mu _{\phi _1}\rangle$ of the order of 3 to 5 mas yr⁻¹. Note that the decrease of the proper motions toward ϕ₁ = 0 coincides with the distance increase of the stream (see Figure 6).

Zoom In Zoom Out Reset image size

Having determined the stream proper motions, we can further improve the CMD-filtered map of the GD-1 stream (Figure 7) by requiring that the proper motions of the stars are consistent with the proper motion of the stream. That is shown in Figure 10 and discussed in Section 3.2.

Zoom In Zoom Out Reset image size

3.2. Radial Velocities

To construct the 6D phase-space distribution of the stream, the radial velocities are the remaining datum. By necessity, the actual sample for which spectra and hence radial velocities will be available will differ from the photometric sample just described. In this section, we will use both the data from the SDSS/SEGUE survey (Yanny et al. 2009) as well as radial velocities obtained by us with the TWIN spectrograph on Calar Alto, specifically targeting likely stream member stars.

3.2.1. SDSS Radial Velocities

SEGUE and SDSS only provide sparse spatial sampling of high latitude stars. SEGUE did not target any GD-1 member stars specifically. Therefore, we have to search through the existing SEGUE spectra to identify likely or possible member stars by position on the sky, CMD position, and proper motion. In the previous section, we described that we used the ratio of the stream/background probabilities $\frac{P_{\rm stream}(u-g,g-r,r-i,r)}{P_{BG}(u-g,g-r,r-i,r)}$ to select high probability members of the stream. Now addition to that, we also select the stars within the μ_α, μ_δ box (see Figure 8). That allows us to have a sample of stream stars with much less background contamination, although the overall size of that sample is significantly smaller, since the SDSS/USNO-B1.0 measurements of the proper motions were done for stars with r ≲ 20 (Munn et al. 2004). To illustrate how good the proper motion selection is when combined with the color–magnitude selection, we show map of high probability stream member stars in Figure 10. The stream is now clearly seen in individual stars. In Figure 10, we also overplot the location of existing SEGUE DR7 pointings, some of which cover the stream. Therefore, we may expect to find some stream members among the SEGUE targets in these fields.

Figure 11 shows the SDSS/SEGUE radial velocity distribution as a function of ϕ₁ for those stars whose proper motions and color–magnitude position are consistent with stream membership, and which are located within 3° from the center of the stream. The typical uncertainty of the SDSS/SEGUE radial velocities is ∼20 km s⁻¹. The filled red circles in this figure show the subset of stars located within 03 from the center of the stream and therefore represent the subset with very high membership probability, while the open black circles represent (spatially selected, |ϕ₂|>03) background stars with similar proper motion and color magnitude. The filled symbols in Figure 11 clearly delineate the radial velocity of the stream. Clumps of red circles are visible at (ϕ₁, V_rad) ≈ (−25°, −100 km s⁻¹), (−30°, −80 km s⁻¹), (−47°, 0 km s⁻¹), (−55°, 40 km s⁻¹). In order to perform the formal measurements of the radial velocities, we performed a maximum likelihood fit by a model consisting of two Gaussians: one (wide) Gaussian representing the background distribution, and a second (narrow) modeling velocity distribution of stream stars. This fit gave us the following results: (ϕ₁, V_rad) = (−56°, 39 ± 14 km s⁻¹), (−47°, −7 ± 10 km s⁻¹), (−28°, −61 ± 6 km s⁻¹), and (−24°, −83 ± 9 km s⁻¹).

Zoom In Zoom Out Reset image size

3.2.2. Calar Alto Radial Velocities

Since the SDSS/SEGUE radial velocities only provide constraints at a few points of the stream, we decided to obtain additional radial velocity information with targeted observations. We based the target selection for likely stream members on all the available information discussed in the previous sections: position on the sky, color–magnitude location, and proper motions. We selected 34 stars likely members with r≲ 19 for the observations. All these stars are within the sample plotted in Figure 10. These stars are mostly main-sequence and main-sequence turn-off stars with 18 ≲ r ≲ 19.

The observations were performed by the TWIN spectrograph on the 3.5 m telescope at Calar Alto observatory, during several nights of service observing in 2009 February. The TWIN spectrograph is the intermediate resolution spectrograph installed in the Cassegrain focus of the telescope. It consists of two separate spectroscopic channels (blue and red) behind the common entrance slit. The light from the slit is split into "red" and "blue" beams by a dichroic mirror. We used the blue and red arms of the spectrograph at a resolution of 4000–5000 to observe the H_β, Mgb lines and Ca ii near-IR triplet, respectively. The standard data reduction steps were applied to the data set using custom routines written in Python language. The median signal-to-noise ratio (S/N) per pixel was 7 for the blue spectra and 3 for red spectra.

We used both the blue and the red spectra to compute the radial velocity of each star. The radial velocity of each star was derived by minimizing χ² as a function of velocity shift of the template convolved with the appropriate Line Spread Function. The χ² for each star was a sum of the χ² for the blue and the red part of the spectra. As the template in the blue spectral range, we used the spectra from the ELODIE database (Prugniel et al. 2007) for stars of similar color and magnitude to the targeted ones and with low metallicity [Fe/H] ∼ −2. In the red spectral range, the template spectrum was simply consisting of three lines of Ca triplet at 8498.02 Å, 8542.09 Å, and 8662.14 Å. The error of each velocity measurement was determined using the condition Δ(χ²(V)) = 1.

The Calar Alto measurements of the velocities together with their errors are overplotted with blue symbols in Figure 11. It is apparent from Figure 11 that for −50° < ϕ₁ < −10° the targeting strategy was very successful, nicely delineating the projected velocity gradient along the stream. Overall out of 34 observed stars, ∼ 24 stars belong to the stream and five did not have enough S/N for the velocity determination. Unfortunately, the targeting near ϕ₁ ≈ +5° failed to identify stream members, probably because the stream there is less intense and further away.

In the previous section, we described the derivation of different stream properties such as distance, position on the sky, and proper motion separately. In this section, we will map the stream in 6D position-velocity space in a more consistent way, using all the available information (see, e.g., Cole et al. 2008, for the application of similar, although simpler method to Sgr stream). This will provide us with a set of orbit constraints along different sections of the GD-1 stream, which we will then model by an orbit to derive potential constraints.

4.1. Positions on the Sky and Distances to the Stream

We start by characterizing the projected stream position and its distance from the Sun through a maximum likelihood estimate for a parameterized model of the stream P_stream(r, g − r, ϕ₁, ϕ₂) that describes it in four-dimensional space of photometric observables r, g − r, ϕ₁, ϕ₂:

$$\begin{eqnarray} &&P_{\rm stream}(r,g-r,\phi _1,\phi _2) = {\rm CMD}(r,g-r,D(\phi _1)) \nonumber \\ &&\quad\times I(\phi _1) \times \frac{1}{\sqrt{2\,\pi }\sigma _{\phi _2}(\phi _1)} {\rm exp}\left(-\frac{(\phi _2-\phi _{2,0}(\phi _1))^2}{2\,\sigma _{\phi _2} ^2(\phi _1) } \right). \end{eqnarray} \tag{ 2 }$$

Here ϕ_2,0(ϕ₁) is the ϕ₂ position of the stream center on the sky as a function of ϕ₁, $\sigma _{\phi _2}(\phi _1)$ is the projected width of the stream in ϕ₂, I(ϕ₁) is the "intensity" (i.e., the number density) of the stream as a function of ϕ₁, and D(ϕ₁) is the distance to the stream. Further, CMD(r, g − r, D(ϕ₁)) is the normalized Hess diagram (i.e., the probability distribution in CMD space) expected for the stream's stellar population at a distance of D(ϕ₁) after accounting for the observational errors. We construct that CMD based on the age and metallicity obtained in Section 2 and the isochrones from Girardi et al. (2000) and Marigo et al. (2008), assuming that $\frac{d[Fe/H]}{d\phi _1}=0$ and $\frac{d(age)}{d\phi _1}=0$. In this model, P_stream depends on four functions—I(ϕ₁), ϕ_2,0(ϕ₁), $\sigma _{\phi _2}(\phi _1)$, and D(ϕ₁)—which we take to be piecewise constant; i.e. for intervals δϕ₁ they simply become four parameters.

For the field stars P_BG(r, g − r, ϕ₁, ϕ₂), the analogous four-dimensional (4D) distribution is

$$$ P_{BG}(r,g-r,\phi _1,\phi _2)= I_{BG}(\phi _1,\phi _2)\times {\rm CMD}_{BG}(r,g-r,\phi _1), $$$

where I_BG(ϕ₁, ϕ₂) is the 2D number density distribution of the field stars around the stream and CMD_BG(r, g − r, ϕ₁) the corresponding CMD. These functions are determined empirically from the data in parts of the sky adjacent to the stream (|ϕ₂| ≳ 05). I_BG(ϕ₁, ϕ₂) is determined by fitting the density of the stars in the ϕ₁, ϕ₂ space by a polynomial. CMD_BG(r, g − r, ϕ₁) is determined by constructing the Hess diagrams using all the stars with 03 < |ϕ₂| < 5° in several ϕ₁ bins.

To simplify the determination of P_BG(r, g − r, ϕ₁, ϕ₂) and P_stream(r, g − r, ϕ₁, ϕ₂), we split the stream in several ϕ₁ pieces, and consider I(ϕ₁), $\sigma _{\phi _2}$, D(ϕ₁), and ϕ_2,0(ϕ₁) as constants within them. The log-likelihood for the mixture of the P_stream and P_BG distribution can be written as (here for convenience, we introduce α as a fraction of stream stars instead of I(ϕ₁)):

$$\begin{eqnarray} {\rm ln}(\mathcal L) &=& \sum \limits _{\rm stars} 2\,{\rm ln} (\alpha \, P_{\rm stream}(r_i,g_i-r_i,\phi _{1,i},\phi _{2,i}) \nonumber \\ &&+(1-\alpha)P_{BG}(r_i,g_i-r_i,\phi _{1,i},\phi _{2,i})) \end{eqnarray} \tag{ 3 }$$

and should be maximized with respect to the parameters $\sigma _{\phi _2}$, D(ϕ₁), ϕ_2,0(ϕ₁) and α. The maximization is performed using the Truncated Newton method (Nash 1984). The parameter errors are obtained using numerically computed Fisher information (see, e.g., Cox & Hinckley 1974).

The left and central panels of Figure 12 show the resulting estimates of the projected position and the distance of the stream, the parameters used in the subsequent orbit fitting. We do not use the number density of stream stars, as it varies noticeably along the stream (see Figure 7) and the reason of these variations is not clear. It is apparent from Figure 10 that the projected stream position is very well defined, and that a distance gradient exists along the stream.

Zoom In Zoom Out Reset image size

4.2. Proper Motions

The likelihood maximization just described also results in stream membership probabilities for any given star i, P_stream(r_i, g_i − r_i, ϕ_1,i, ϕ_2,i). This information can then be used to estimate via maximum likelihood the mean proper motions of different stream pieces, thereby extending the observational estimates to the full 6D space (r, g − r, ϕ₁, ϕ₂, $\mu _{\phi _1}$, $\mu _{\phi _2}$).

$$\begin{eqnarray} &&P_{{\rm stream},\mu }(r,g-r,\phi _{1},\phi _{2},\mu _{\phi _1},\mu _{\phi _2}) \nonumber \\ &&\quad=P_{\rm stream}(r,g-r,\phi _1,\phi _2) \times \frac{1}{2\,\pi \,\sigma ^2_\mu }\nonumber \\ &&\quad\times{\rm exp}\left(-\frac{(\mu _{\phi _1}-\mu _{\phi _1,0}(\phi _1))^2+(\mu _{\phi _2}-\mu _{\phi _2, 0 } (\phi _1))^2 }{ 2\,\sigma ^2_\mu } \right), \end{eqnarray} \tag{ 4 }$$

where we simply take the previously determined P_stream as a prior, to be modified by the Gaussian distribution in the 2D proper motions space. Here, the proper motion distribution is characterized by three functions, $\mu _{\phi _1,0}(\phi _1)$, $\mu _{\phi _2, 0 } (\phi _1),$ and σ_μ, which again we take to be piecewise constant. The distribution of the background stars in the 6D space $P_{BG,\mu }(r,g-r,\phi _1,\phi _2,\mu _{\phi _1},\mu _{\phi _2})$ is obtained empirically by binning the observational data. Then we construct again the logarithm of likelihood, considering variations in four parameters: the number of stream stars, the proper motion of the stream in ϕ₁ and ϕ₂, and the proper motion spread σ_μ. This likelihood is then maximized and we determine the $\mu _{\phi _1}$, $\mu _{\phi _2},$ and σ_μ for different stream pieces. The right panel of Figure 12 shows the resulting proper motion estimates as a function of ϕ₁. These proper motions have not been corrected here for the Sun's reflex motion, which we will model in Section 5.

Overall, the analysis presented in the previous sections has resulted in the best and most extensive set of 6D phase-space coordinate map for a cold stream of stars in the MW.

If we can assume that all the stream stars lie close to one single test-particle orbit, then our phase-space map of the GD-1 stream should not only define this orbit, but at the same time constrain the Milky Way's potential. In general, the stars in the tidal stream have slightly different energies and angular momenta, but the assumption that the stream stars are moving along the same orbit is plausible, especially if the stream is thin and is near pericenter (Dehnen et al. 2004; W. Dehnen 2004, private communication). But it is not straightforward to quantify the quality of such an approximation. For now, we simply fit an orbit to our 6D map of available observational data: the position on the sky, ϕ₂(ϕ₁), proper motion $\vec{\mu }(\phi _1)$, distance to the stream D(ϕ₁), and radial velocity V_rad(ϕ₁). For each assumed potential, we will determine the best fit orbit, but then marginalize over the orbits to determine the range of viable gravitational potentials. This analysis extends earlier efforts by Grillmair & Dionatos (2006b) and Willett et al. (2009) who have presented orbit solutions for GD-1. However, we can now draw on a much more extensive set of observational constraints. We also explore the fit degeneracies. Given that our 6D phase-space map of the GD-1 stream spans only a limited range in R and z (as seen from the Galactic center 11 kpc ≲ R ≲ 14 kpc, 5 kpc ≲ z ≲ 9 kpc), it proved useful to consider very simple parameterized potentials at first. Further, it proved necessary to consider what prior information we have on the Sun's (i.e., the observers) position and motion, as well as on our MW's stellar disk mass.

5.1. One Component Potential

The stream is located at Galactocentric (R, z) ≈ (12, 6) kpc, a regime where presumably both the stellar disk and the dark halo contribute to the potential and its flattening. Of course, the stream dynamics are solely determined by the total potential, and therefore we consider first a simple single-component potential, the flattened logarithmic potential

$$\begin{equation} \Phi (x,y,z)=\frac{V_c^2}{2}\, {\rm ln}\left(x^2+y^2+\left(\frac{z}{q_{\Phi }}\right)^2\right)\mbox{,} \end{equation} \tag{ 5 }$$

which has only two parameters: the circular velocity V_c and the flattening q_Φ. Note that (1 − q_density) ≈ 3(1 − q_Φ) for moderate flattening (e.g., p. 48 of Binney & Tremaine 1987). Such a simple potential seems justified as the stream stars are only probing a relatively small range in R and z.

In practice, we fit an orbit to the 6D stream map by considering a set of trial starting points in the Galaxy, specified by the initial conditions $(\vec{X}(0),\dot{\vec{X}}(0))$ in standard Cartesian Galactic coordinates. Together with an assumed gravitational potential, this predicts $(\vec{X}(t),\dot{\vec{X}}(t))$, which can be converted to the observables, ϕ₂(ϕ₁), $\vec{\mu }(\phi _1)$, D(ϕ₁), V_rad(ϕ₁) and then compared to the 6D observations (Figure 13). For each $[(\vec{X}(0),\dot{\vec{X}}(0)|\Phi (\vec{X})]$ we can evaluate the quality of the fit by calculating χ², summing over all data points, shown in Figure 13. For any given $\Phi (\vec{X})$, χ² can be then minimized with respect to the orbit, i.e., $(\vec{X}(0),\dot{\vec{X}}(0))$, providing the "best-fit" orbit in this potential and the plausibility of that potential. The minimization is performed using the MPFIT code (Markwardt 2009) implementing the Levenberg–Marquardt technique (Marquardt 1963) translated into Python.⁶ The data used to constrain the potential are given in the Tables 1, 2, 3, and 4 (except the SDSS measurements of the radial velocities which are given in the end of Section 3).

Zoom In Zoom Out Reset image size

Table 1. Radial Velocities from the Calar Alto Observations

Star	ϕ1	ϕ2	Vrad
	(deg)	(deg)	(km s−1)
SDSS J094105.35+315111.6	−45.23	−0.04	28.8 ± 6.9
SDSS J094705.26+332939.8	−43.17	−0.09	29.3 ± 10.2
SDSS J095740.48+362333.0	−39.54	−0.07	2.9 ± 8.7
SDSS J095910.43+363206.6	−39.25	−0.22	−5.2 ± 6.5
SDSS J100222.01+374113.3	−37.95	0.00	1.1 ± 5.6
SDSS J100222.02+374049.2	−37.96	−0.00	−11.7 ± 11.2
SDSS J101033.02+393300.8	−35.49	−0.05	−50.4 ± 5.2
SDSS J101110.08+394453.9	−35.27	−0.02	−30.9 ± 12.8
SDSS J101254.83+395525.6	−34.92	−0.15	−35.3 ± 7.5
SDSS J101312.05+400613.3	−34.74	−0.08	−30.9 ± 9.2
SDSS J101702.15+404747.3	−33.74	−0.18	−74.3 ± 9.8
SDSS J101951.76+412701.5	−32.90	−0.15	−71.5 ± 9.6
SDSS J102216.20+415534.7	−32.25	−0.17	−71.5 ± 9.2
SDSS J103003.87+434351.7	−29.95	−0.00	−92.7 ± 8.7
SDSS J104341.92+460224.7	−26.61	−0.11	−114.2 ± 7.3
SDSS J104840.98+464922.1	−25.45	−0.14	−67.8 ± 7.1
SDSS J105036.96+472000.1	−24.86	0.01	−111.2 ± 17.8
SDSS J110711.27+494415.9	−21.21	−0.02	−144.4 ± 10.5
SDSS J114242.08+533841.4	−14.47	−0.15	−179.0 ± 10.0
SDSS J114724.59+535546.8	−13.73	−0.28	−191.4 ± 7.5
SDSS J115116.08+542142.7	−13.02	−0.21	−162.9 ± 9.6
SDSS J115326.06+542930.6	−12.68	−0.26	−217.2 ± 10.7
SDSS J115404.06+543511.4	−12.55	−0.23	−172.2 ± 6.6

Download table as:  ASCIITypeset image

Table 2. Stream Positions

ϕ1	ϕ2
(deg)	(deg)
−60.00	−0.64 ± 0.15
−56.00	−0.89 ± 0.27
−54.00	−0.45 ± 0.15
−48.00	−0.08 ± 0.13
−44.00	0.01 ± 0.14
−40.00	−0.00 ± 0.09
−36.00	0.04 ± 0.10
−34.00	0.06 ± 0.13
−32.00	0.04 ± 0.06
−30.00	0.08 ± 0.10
−28.00	0.03 ± 0.12
−24.00	0.06 ± 0.05
−22.00	0.06 ± 0.13
−18.00	−0.05 ± 0.11
−12.00	−0.29 ± 0.16
−2.00	−0.87 ± 0.07

Download table as:  ASCIITypeset image

Table 3. Stream Distances

ϕ1	Distance
(deg)	(kpc)
−55.00	7.20 ± 0.30
−45.00	7.59 ± 0.40
−35.00	7.83 ± 0.30
−25.00	8.69 ± 0.40
−15.00	8.91 ± 0.40
0.00	9.86 ± 0.50

Download table as:  ASCIITypeset image

Table 4. Stream Proper Motions

ϕ1	$\mu _{\phi _1}$	$\mu _{\phi _2}$	σμ
(deg)	(mas yr−1)	(mas yr−1)	(mas yr−1)
−55.00	−13.60	−5.70	1.30
−45.00	−13.10	−3.30	0.70
−35.00	−12.20	−3.10	1.00
−25.00	−12.60	−2.70	1.40
−15.00	−10.80	−2.80	1.00

Download table as:  ASCIITypeset image

It is crucial to note that the conversion of $(\vec{X}(t),\dot{\vec{X}}(t))$ to the space of observables depends on the position and motion of the observer, i.e., on distance from the Sun to the Galactic center (R₀) and on the 3D velocity of the Sun in the Galaxy rest frame ($\vec{V}_0$). At this stage, we adopt R₀ = 8.5 kpc based on recent determinations (e.g., Ghez et al. 2008; Gillessen et al. 2009), but later we will relax this. The second parameter $\vec{V}_0\equiv \vec{V}_{{\rm LSR}} + \Delta \vec{V}_{{\rm LSR}}$ (where V_LSR is the velocity of the local standard of rest (LSR) and $\Delta \vec{V}_{{\rm LSR}}$ is the Sun's velocity relative to the LSR) is linked to the fitting not only through conversion of the observable relative stream velocities to the GC rest system, but also conceptually through the plausible demand that $\Phi (\vec{X})$ and in particular V_c(R₀, 0) also reproduces $\vec{V}_{{\rm LSR}}$. In this way, constraints on the potential flattening can be derived by considering $r\frac{d\Phi }{dr}$ in the disk plane ($\vec{V}_{{\rm LSR}}$) and the plane of the GD-1 stream. The velocity of the Sun relative to the LSR $\Delta \vec{V}_{{\rm LSR}}$ is quite well known from the HIPPARCOS measurements (Dehnen & Binney 1998): $\Delta \vec{V}_{{\rm LSR}}[{\rm km \;s^{-1}}]=10\vec{e}_x+5.25\vec{e}_y+7.17\vec{e}_z$. The velocity of the LSR, i.e., V_c(R₀, 0) has a considerable uncertainty (Brand & Blitz 1993; Ghez et al. 2008; Xue et al. 2008; Reid et al. 2009). Initially, we will consider the velocity of the LSR simply a consequence of the assumed potential, i.e., V_LSR ≡ V_c(R₀, 0).

Figure 13 illustrates the result of such fitting, by overplotting the best-fit orbit for the plausible potential with V_c = 220 km s⁻¹ and q_Φ = 0.9 over observational data. It is clear that even for the simple flattened logarithmic potential, an orbit can be found that reproduces the observables well. This fit serves to illustrate a few generic points that also hold for orbit fits in differing potentials: the stream moves on a retrograde orbit and it is near pericenter, where the stream is expected to approximate an orbit well (Dehnen et al. 2004). After fitting a first orbit, we may also note its global parameters (see Figure 14 for a 3D map of the orbit): the pericenter is at 14 kpc from the GC; the apocenter is at 26 kpc; and the inclination is 39°.

Zoom In Zoom Out Reset image size

For any given potential $\Phi (\vec{X}|V_c, q_\Phi)$ we can find best-fit orbit by marginalizing over $(\vec{X}(0),\dot{\vec{X}}(0))$ to see what our 6D map of GD-1 implies about the relative plausibility of different V_c and q_Φ. Figure 15 shows the log-likelihood surface for the potential parameters (V_c, q_Φ); note again that this fit neglects all other prior information on V_c at the Sun's position. The contours show 1σ, 2σ, and 3σ confidence regions on the parameters, derived from the $\delta ({\rm ln}(\mathcal {L}))$ values for 2 degrees of freedom (i.e., a two parameter fit; see, e.g., Lampton et al. 1976). The insets at the left and bottom show the marginalized distributions for single parameters. The best-fit values with the two-sided 68% confidence intervals are V_c = 221⁺¹⁶₋₂₀ km s⁻¹ and q_Φ = 0.87^+0.12_−0.03. Figure 15 illustrates that the flattening parameter q_Φ is quite covariant with the equatorial circular velocity V_c. An extreme example may serve to explain this covariance qualitatively. If the stream went right over the pole (z-axis), then the local force gradient would be proportional to V_c × q_Φ (Equation (5)). Information about the potential flattening must therefore come from combining kinematics and dynamics in the disk plane with the information from GD-1.

Zoom In Zoom Out Reset image size

The fit shown in Figure 15 asks the data not only to constrain the potential at the stream location and determine the stream orbit, but also to infer the Sun's motion (or at least V_LSR) from its reflex effect on the data. Clearly, providing a prior on V_c(R₀, 0) is sensible, especially if we care about constraints on the potential flattening. We consider the constraints that arise from the Sun's reflex motion with respect to the Galactic center, the most robust and sensible prior in this context. Ghez et al. (2008) recently combined radio data (Reid & Brunthaler 2004) with near-IR data on the Galactic center kinematics to arrive at V_c(R₀, 0) = 229 ± 18 km s⁻¹. It is also noticeable that our own constraint on V_c(R₀, 0) from Figure 15, 221⁺¹⁶₋₂₀ km s⁻¹, is close both in value and uncertainty to the estimate of Ghez et al. (2008), which is based on a completely disjointed data set and approach.

Figure 16 shows the resulting log-likelihood contours and 1σ, 2σ, and 3σ confidence regions after applying this prior on the V_c (229 ± 18 km s⁻¹). Note that likelihood on Figure 16 is also marginalized over R₀ with a Gaussian prior (R₀ = 8.4 ± 0.4 kpc Ghez et al. 2008).

Zoom In Zoom Out Reset image size

Figure 16 shows that the posterior probability distribution on V_c has slightly changed to V_c(R₀) = 224⁺¹²₋₁₄ km s⁻¹ with noticeably smaller error bar compared to the value from Ghez et al. (2008). Figure 16 also shows us the slight improvement comparing to Figure 15 of flattening constraints: q_Φ = 0.87^+0.07_−0.04. This means that the total potential appears to be oblate (in the radial range probed); this may not be surprising, as the stellar disk—which is manifestly very flattened—contributes to the total potential.

Figure 17 illustrates how well the best-fit orbits for different potentials $\Phi (\vec{X}|V_c,q_\Phi)$ can mimic one other in the space of observables. This is the source of the parameter covariances shown in Figures 15 and 16.

Zoom In Zoom Out Reset image size

The fitting of the orbit shown in Figure 13 allows us to make an estimate of the line-of-sight velocity dispersion in the stream, by comparing the dispersion of the radial velocity residuals with the accuracy of individual radial velocities. Assuming that without observational uncertainties the residuals from the orbit fit should be Gaussian distributed with zero mean, we performed the maximum likelihood fitting of the residuals, which gives a 90% confidence upper limit on the velocity dispersion of stars in the stream of ∼3 km s⁻¹.

Before we aim at separating possible flattening contributions from the halo and disk, it is worth commenting on the accuracy and limitations of our estimate of q_Φ. In the range (R, z) ≈ (12, 6) kpc no other direct observational constraints on the potential shape exist in the literature, and hence our estimate of q_Φ = 0.87^+0.07_−0.04 is a new and important contribution. On the other hand, an error of δq_Φ ∼ 0.05, especially when translated into the flattening error of the equivalent scale-free mass distribution, may not appear as particularly helpful in model discrimination, or as impressively accurate—especially as a manifestly cold stellar stream spanning over 60° on the sky may seem ideal for mapping the potential at first glance.

5.2. Constraints on the Shape of the Dark Matter Halo from a Bulge, Disk, Halo Three-component Potential

In the previous section, we constrained the parameters of a simplified MW potential, the spheroidal logarithmic potential. It is clear that the MW potential at the position of the stream must depend explicitly on the sum of baryonic Galaxy components (bulge and disk) and on the dark matter halo. We now explore whether our constraint on the shape of the overall potential, q_Φ ∼ 0.9, permits interesting statements about the shape of the DM potential itself. At the distance of (R, z) ≈ (12, 6) kpc, the contribution of the disk to the potential should still be relatively large, weakening or at least complicating inferences on the shape of the DM distribution.

We adopt a three-component model of the Galaxy potential, choosing one that is widely used in the modeling of the Sgr stream (Helmi 2004; Law et al. 2005; Fellhauer et al. 2006) and reproduces the galactic rotation curve reasonably well.

The model consists of a halo, represented by the logarithmic potential The numerical values of parameters of the potential were taken from Fellhauer et al. (2006)

$$\begin{equation} \Phi _{\rm halo}(x,y,z)=\frac{v_{\rm halo}^2}{2}\, {\rm ln}\left(x^2+y^2+\left(\frac{z}{q_{\Phi,\rm halo}}\right)^2+d^2\right)\mbox{,} \end{equation} \tag{ 6 }$$

where we have adopted d = 12 kpc from the previous authors. The disk is represented by a Miyamoto–Nagai potential (Miyamoto & Nagai 1975),

$$\begin{equation} \Phi _{\rm disk}(x,y,z)=\frac{G M_{\rm disk}}{\sqrt{x^2+y^2+(b+\sqrt{z^2+c^2})^2}} \end{equation} \tag{ 7 }$$

with b = 6.5 kpc, c = 0.26 kpc. The bulge is modeled as a Hernquist potential:

$$\begin{equation} \Phi _{\rm bulge}(x,y,z)=\frac{G M_b}{r+a} \end{equation} \tag{ 8 }$$

with M_b = 3.4 × 10¹⁰ M_☉, a = 0.7 kpc.

As in the previous section, for any given set of potential parameters, we can find the best-fitting stream orbit and compute χ² of the fit. In the current paper, we do not make any attempts to fully fit all the parameters of the MW potential, but we try to make an estimate of the MW DM halo flattening. We take the three-component potential and fix all but three parameters—disk mass M_disk, circular velocity of the halo v_halo, and the flattening of the halo q_Φ,halo. On a 3D grid of these parameters, we perform a χ² fit. Figure 18 shows the results of such fit after marginalization over the orbital parameters $(\vec{X}(0),\dot{\vec{X}}(0))$, circular velocity of the halo v_halo with a Gaussian prior from Xue et al. (2008) and a Gaussian prior on the circular velocity at the Sun's radius from Ghez et al. (2008). The figure clearly illustrates that in the case of the three-component potential, the current data are unable to give a significant new insights on the flattening of the DM halo. We can only say that at 90% confidence q_Φ,halo > 0.89. We note however that for the future analysis, if a multi-component model for the potential is used then more prior information is required, i.e., not only on V_LSR and v_halo but also on M_disk and other parameters of the potential.

Zoom In Zoom Out Reset image size

In this paper, we have presented a thorough analysis of the GD-1 stream combining the publicly available SDSS and SEGUE data with follow-up spectroscopic observations from Calar Alto. The combination of the photometric SDSS observations, USNO-B/SDSS proper motions, SDSS, SEGUE, and Calar Alto radial velocities allowed us to construct a unprecedented 6D phase-space map of this very faint (29 mag/sq. arcsec) tidal stream. The 6D phase-space map of the stream, spanning more than 60° on the sky, provided the opportunity not only to fit the orbit as Grillmair & Dionatos (2006b) and Willett et al. (2009) have done previously but also to explore what constraints can be placed on the MW potential.

The analysis is based on the assumption that the stream stars occupy one orbit. In detail, of course, different stars on the stream have slightly different values of conserved quantities and therefore lie on slightly different orbits. Effectively, our analysis depends on these being small when compared with an orbital uncertainties. The magnitude of the departure of the stream from a single orbit will, in detail, be a function of the progenitor and the disruption process; as these details became understood the model can be refined.

We found that the distribution of stream stars in phase space can be well fit by an inclined eccentric orbit in the spheroidal logarithmic potential. After marginalization over the stream orbital parameters, we derive a circular velocity V_c = 221⁺¹⁶₋₂₀ km s⁻¹ and flattening q_Φ = 0.87^+0.12_−0.03. This measurement has been made without the use of any information other than that in the GD-1 stream itself. It is important that the information available in the observations of the stream is very sensitive to the V_c, the circular velocity at the Sun's position. The reason for that is that the stream extends more then 60° on the sky and therefore both the radial velocities and the proper motions have components coming from the projection of the Sun's motion.

If we combine our circular velocity measurement with existing prior on the V_c from Ghez et al. (2008) and also marginalize over the distance from the Sun to the Galactic center using the Ghez et al. (2008) prior (R₀ = 8.4 ± 0.4 kpc), we further tighten the error bar on V_c = 224⁺¹²₋₁₃ km s⁻¹ and on the flattening of the potential q_Φ = 0.87^+0.07_−0.04. Our measurement of the V_c is the best constraint to date on the circular velocity at the Sun's position, and the measurement of q_Φ is the only strong constraint on q_Φ at galactocentric radii near R ∼ 15 kpc.

The measurement of the flattening of the potential q_Φ = 0.87^+0.07_−0.04 describes only the flattening of the overall Galaxy potential at the stream's position (R, z) ≈ (12, 6) kpc where the disk contribution to the potential is presumably large. Unfortunately the data on the GD-1 stream combined with the Ghez et al. (2008) and Xue et al. (2008) priors on V_c and v_halo are not enough for separating the flattening of the halo from the flattening of the total Galaxy potential. So, we are unable to place strong constraints on flattening of the MW DM halo; we put a 90% confidence lower limit at q_Φ,halo > 0.89.

Despite the negative result on the measurement of the MW DM halo flattening, we note that the data from the GD-1 stream is able to give strong constraints on two important Galaxy parameters. We claim that that our data set on the GD-1 stream should now be combined with other available MW kinematical data (from other stellar streams, blue horizontal branch (BHB) stars, MW rotation curve) in order to tighten the existing constraints on the Galaxy parameters. It is important that the constraints on Galaxy potential based on the GD-1 stream data set are, to large extent, model-independent and purely kinematic, i.e., the constraints on the V_c come to large extent from the projection effects and manifest themselves in proper motions and radial velocities.

We think that it is quite surprising that such a long (60°) stream with full 6D map did not allow us to constrain large number of parameters of the MW DM halo and other Galactic components. We think that there are several reasons for that, but the most important is that while the observed part of the orbit spans ∼ 70° on the sky from the Sun's point of view, the orbital phase spanned by the stars as seen from the Galactic center is only ∼40°. Furthermore, the observable part of the stream occupies the perigalacticon part of the orbit, so the range of Galactocentric distances probed by the stream is small. This makes it plausible that orbits of different eccentricities, and hence of different azimuthal velocities at their pericenter, can match so closely the same set of 6D coordinates. Since very cold streams take many orbits to spread a substantial fraction of 2π in orbital phase (see, e.g., also Pal 5; Odenkirchen et al. 2001; Grillmair & Dionatos 2006a), all future analyses of yet-to-be discovered streams will particularly need to consider the trade-off between the conceptual and practical attractiveness of "cold" streams and the near-inevitable limitations of their phase coverage.

In addition to the weakness coming from phase coverage, our analysis at this point must rely on photometric distance estimates; these have random errors of ∼10%, after an empirical distance correction to the best-fit isochrones that is of the same magnitude (see Section 2). The proper motions that we derived for ensembles of stream stars are unprecedented for a stellar stream in the MW's outer halo; however, the corresponding velocity precision, especially when compounded by distance errors, is still the largest single source of uncertainty in the fitting (e.g., our tests have shown that overestimating heliocentric distances to the stream leads to the overestimated measurement of V_c). Both the deficiencies of the distances and proper motions on GD-1 will be largely alleviated after the launch of the GAIA satellite, which will allow much tighter constraints on the Galactic potential.

Before the launch of GAIA, we suggest that further observations of radial velocities of stars in the stream (with deeper spectroscopy) and improvements in proper motion precision (with, e.g., Pan-STARRS; Kaiser et al. 2002) should be able to reduce the error bars on Galaxy parameters significantly. It is also important to properly calibrate the distance to the stream, which may be done by confirming several probable BHB candidates in the stream.

We suggest that any attempt to fit the Galactic potential (such as Widrow et al. 2008) now should not ignore the data set on GD-1 and should incorporate it into their fits.

The orbital parameters, which we have measured for the GD-1 stream, are more or less consistent with those from Willett et al. (2009): pericenter is at 14 kpc from the GC, the apocenter is at 26 kpc, and the orbit inclination with respect to the Galactic plane is 39°. We have also estimated the total stellar mass associated with the stream to be ∼2 × 10⁴ M_☉, which together with the relatively small stream width of ∼20 pc suggests that the progenitor of the stream was a globular cluster, although we cannot completely rule out the dwarf galaxy progenitor. We have also determined the 90% confidence upper limit for the velocity dispersion of the stars in the stream to be ∼3 km s⁻¹, which does not significantly contradict either the globular cluster or dwarf galaxy progenitor hypothesis. Given the length of the observable part of the stream of ∼10 kpc and the velocity dispersion ≲3 km s⁻¹, the age of the stream can be estimated to be greater than 1.5 Gyr (assuming that the progenitor is in the middle of the observed part of the stream). It is interesting that the stream managed to evade possible destruction by interaction with DM subhalos orbiting around MW (Carlberg 2009). Although, the clumpiness observed in the stream may be attributed to these past interactions (S. E. Koposov et al. 2010, in preparation).

Overall in this paper, we have illustrated a method for analyzing the thin stellar stream using all the available information on it and further utilizing that information to constrain the Galaxy potential. We believe that in the epoch of Pan-STARRS, LSST and especially GAIA, which will give us a wealth of new information on the MW halo, stream-fitting like that presented in this paper will be extremely useful and productive.

S.K. was supported by the DFG through SFB 439, by a EARA-EST Marie Curie Visiting fellowship, and partially by RFBR 08-02-00381-a grant. S.K. acknowledges hospitality from the Kavli Institute for Theoretical Physics (KITP) Santa Barbara during the workshop "Building the Milky Way." S.K. thanks Jelte de Jong for running MATCH code on the GD-1 data and the Calar Alto observing staff for excellent support. D.W.H. acknowledges support from NASA (grant NNX08AJ48G) and a Research Fellowship from the Alexander von Humboldt Foundation. The authors thank Kathryn Johnston, Jorge Peñarubia, and James Binney for useful discussions and the anonymous referee for the elaborated referee report which helped us improve the paper.

The project made use the open-source Python modules matplotlib, numpy, scipy; SAI Catalogue Access Services (Sternberg Astronomical Institute, Moscow, Russia; Koposov et al. (2007a); NASA's Astrophysics Data System Bibliographic Services).

The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.