Fitting the Density Substructure of the Stellar Halo with MilkyWay@home

The distribution of stars in the Galactic halo is dominated by dwarf galaxies and tidal streams of stars that have been stripped from them (see Figure 1 of Newberg et al. (2002), and the "Field of Streams" from Belokurov et al. (2006)). These substructures represent the recent minor merger history of the Milky Way (Bullock & Johnston 2005), and they contribute to the buildup of stars in the Milky Way stellar halo. Two or three dozen tidal debris streams, most of which extend tens of degrees or more across the sky, have been identified; the exact number cannot be determined due to controversy over the identity of tidal streams, particularly those discovered near the Galactic plane, and because some streams are detected at low enough significance that they are considered stream "candidates." In addition to streams, other substructures of ambiguous origin (most notably "clouds") continue to be discovered in the Galactic spheroid.

For a review of tidal streams and clouds, see Grillmair & Carlin (2016). A more recent list of streams and clouds included in the "GALSTREAMS" Python Package can be found in Table 4 of Mateu et al. (2018). More recent halo substructures are identified in Li et al. (2016b), Sohn et al. (2016), Grillmair (2017a, 2017b), Jethwa et al. (2017), and Shipp et al. (2018). In particular, Shipp et al. (2018) identify eleven new substructures in the Milky Way halo using data from the Dark Energy Survey (DES; The Dark Energy Survey Collaboration 2005, 2016).

Identification and measurement of tidal debris in the Milky Way halo is useful for understanding structure formation and galaxy assembly, and it has the potential to constrain the density distribution of the Milky Way's stellar halo. Methods for measuring the halo shape from tidal streams have been, and continue to be, developed (e.g., Law & Majewski 2010; Koposov et al. 2013; Küpper et al. 2015; Bovy et al. 2016; Johnston & Carlberg 2016; Dierickx & Loeb 2017a; Sanderson et al. 2017). In addition, the distribution of dark subhalos can be measured by looking for stars ejected from tidal streams (Siegal-Gaskins & Valluri 2008), stream heating (Johnston 2002) or stream gaps (Carlberg 2012). Pearson et al. (2017) show that streams can also be used to constrain the rotation rate of the Galactic bar.

These techniques to determine the dark matter distribution in the Milky Way from tidal streams rely on accurate measurements of the tidal debris itself, but since we discovered that halo tidal streams are more numerous and complex than originally thought, the association of particular stars with particular tidal streams has become more ambiguous. For example, Newberg et al. (2009) discovered that the blue horizontal branch stars (BHBs) thought to be associated with the southern portion of the Sagittarius (Sgr) dwarf tidal stream in Yanny et al. (2000) are actually part of the Cetus Polar Stream. The Sgr dwarf tidal stream, which is the most prominent tidal stream in the sky, and the so-called "bifurcated" Sgr stream that appears to split off from it, have also caused confusion; for example, Newby et al. (2013) suggested that the southern Sgr stream could be associated with the "bifurcated" stream in the north, and the northern Sgr stream could be associated with the "bifurcated" stream in the south. These misidentifications and possible misidentifications of stars in the most prominent halo streams underscore the difficulties in counting and characterizing tidal streams.

In this paper, we present an improved statistical photometric parallax (Cole et al. 2008; Newberg 2013) method to measure the spatial density of stars in the Milky Way stellar halo, using turnoff stars from the Sloan Digital Sky Survey (SDSS; York et al. 2000). Statistical photometric parallax is the use of statistical knowledge of the distribution of the absolute magnitudes of stellar populations to determine the underlying density distributions of those stars. This differs from photometric parallax in that the distance to each individual star is not determined.

The idea of using turnoff stars to trace Milky Way halo substructure was introduced by Newberg et al. (2002). They observed density substructure in the SDSS turnoff stars on the Celestial equator and fit an absolute magnitude distribution to their blue turnoff star tracers. In Cole et al. (2008), these tracers and a simplified absolute magnitude distribution from Newberg et al. (2002) were used to build a model of the Milky Way halo and its substructure, and to complete preliminary fits to the stellar density of the halo and the Sagittarius dwarf galaxy tidal stream. Several years later, Newby et al. (2013) continued working with this model and showed that a massive distributed computing network, MilkyWay@home, could be effective in constraining the parameters in the density models of tidal streams. Taking advantage of this new computational power, several optimizations were run on each stripe. Initially, the fitting algorithms were allowed great freedom in selecting the parameters. Later, the parameters were constrained based on the results from neighboring stripes.

These previous studies successfully used statistical photometric parallax to study the structure of the halo using SDSS turnoff stars. SDSS turnoff stars, detected to a limiting magnitude of g = 22.5, can be used to trace the structure of the Milky Way to 45 kpc from the Sun. However, the turnoff stars in a single stellar population, with the same color, can differ in absolute magnitude by two magnitudes (producing a distance error of a factor of 2.5). Photometric parallax (e.g., Jurić et al. 2008) is unusable with turnoff stars because astronomers do not have a way to determine the distance to individual stars with reasonable accuracy using photometry alone.

It has been shown that the absolute magnitude distribution of turnoff stars in halo globular clusters are surprisingly similar to each other, over a metallicity range −2.3 < [Fe/H] < −1.2 dex and over ages ranging from 9 to 13.5 Gyr (Newby et al. 2011). Grabowski et al. (2013) showed that this similarity holds even for the globular cluster Whiting 1, which is only 6 Gyrs old and has a metallicity of approximately [Fe/H] ∼ 0.6 dex (Carraro et al. 2007; Valcheva et al. 2015). This surprising result, which comes about due to the age–metallicity relation for Milky Way stars, makes turnoff stars very useful for tracing the density of the stellar spheroid and outer disk.

In our work, we improve on the statistical photometric parallax methods by implementing a better model for the absolute magnitude of the tracer stars and their detection efficiency, and by using better fitting methods on MilkyWay@home. In our implementation of statistical photometric parallax, we find the parameters in a density model that make the apparent magnitudes and angular positions of the observed stars most likely, using a maximum likelihood estimator (MLE) (Ivezić et al. 2014). The statistical description of the absolute magnitudes of the stellar tracers, and of the selection effects in the data, make statistical photometric parallax somewhat complex to apply. Because we are able to take all of these effects into account, we can reliably measure density distributions in real data.

There are four parts to statistical photometric parallax: data, a density model, an algorithm for measuring how well the model fits the data, and an algorithm for optimizing parameters. The algorithm that measures how well the model fits the data includes the main sequence turnoff (MSTO) absolute magnitude distribution, as well as any observational biases. It has taken us many years to perfect the algorithm that can simultaneously fit the spatial density of several tidal streams plus a smooth distribution to the SDSS MSTO stars; the smooth distribution represents the sum of streams from small satellites, old streams that are well-mixed in density, and stars that were created during the collapse of the Milky Way, if any.

In this paper, we describe an improved algorithm for characterizing the spatial characteristics of stellar streams in the Milky Way halo using turnoff stars, and show that it is capable of simultaneously recovering the characteristics from three tidal streams plus a smooth halo component, using simulated data designed to mimic the stellar density in the actual Milky Way halo.

We will present preliminary results for one 25 wide SDSS stripe (stripe 19) of data that cuts across the northern Galactic hemisphere. Figure 1 shows the position of stripe 19 in the SDSS northern footprint. The results from this stripe provide measurements of the largest known substructures in the Milky Way halo: the Sgr dwarf tidal stream, the so-called "bifurcated" stream, and the Virgo Overdensity. Stripe 19 crosses the Sgr dwarf tidal stream and the bifurcated stream, in a region of the sky in which they are clearly separated. Given that stripe 19 is more than 20° from the densest portion of the Virgo Overdensity, it is uncertain whether a third substructure measured here is in the tails of the Virgo Overdensity, or whether it is associated with a new halo substructure. In addition to three substructures, we fit smooth Milky Way halo and thick disk distributions.

Zoom In Zoom Out Reset image size

The Sgr dwarf galaxy was first discovered by Ibata et al. (1995), who found evidence of a dwarf galaxy within 16 kpc of the Galactic center, on the far side of the Milky Way, that was thought to be in the process of tidally disrupting. The tidal stream of stars stripped from this dwarf galaxy have since been found to dominate the substructure of the Galactic halo (e.g., Newberg et al. 2002; Majewski et al. 2003; Belokurov et al. 2006; Hernitschek et al. 2017). Though the Sgr dwarf galaxy and the stream of stars that have been tidally stripped from its gravitational grasp have been studied extensively (see Law & Majewski (2016), for a recent review); we are only starting to understand the dynamical history of this present-day merger.

It has been a challenge to find a disruption model that simultaneously fits the positions of the leading and trailing tidal streams in the sky, the line-of-sight velocities of the stream stars, and the observed extension of the trailing tidal tail to ∼100 kpc from the Galactic center (Newberg et al. 2003; Belokurov et al. 2014). Dierickx & Loeb (2017b) present a recent simulation of the tidal debris that reproduces most of the measurements of the position of the leading and trailing tidal debris, including the distant stars in the trailing tidal tail and the observed "spurs" at apogalacticon (Sesar et al. 2017), but it still does not reproduce the line-of-sight velocities of the leading tail. A previous model by Law & Majewski (2010) was able to fit the velocities of the leading tidal tail, using a triaxial dark halo model in which the disks rotate around the intermediate axis. However, this Milky Way configuration is very unlikely (Debattista et al. 2013). Refining the spatial distribution of the Sgr dwarf tidal stream using the algorithm described in this paper will help constrain N-body simulations of the Sgr dwarf tidal disruption and lead to a better understanding of the shape of the Milky Way's dark matter halo.

The "bifurcated" stream can be seen clearly in the "Field of Streams" as a lower surface brightness companion stream to the Sgr dwarf tidal stream (Belokurov et al. 2006). Belokurov identifies the Sgr dwarf tidal stream as Stream A, the "bifurcated stream" as Stream B, and tentatively identifies a more distant Stream C behind Stream A, which is now generally associated with an extension of the Sgr trailing tidal tail (Li et al. 2016a). Koposov et al. (2012) shows the analogous bifurcated stream in the southern Galactic cap. Although the origin of the second, lower surface brightness stream close to the Sgr stream is not known, a leading possibility is that it could arise from multiple wraps of the stream around the Milky Way (Fellhauer et al. 2006). Since its discovery, this stream has remained relatively unstudied compared to its sibling. Newberg et al. (2007) derive distances that are slightly farther away than Sgr for the bifurcated stream. In contrast, Niederste-Ostholt et al. (2010) say that the bifurcated stream is slightly closer to the Sun than the Sgr dwarf tidal stream, and Ruhland et al. (2011) finds the distances to be basically the same. Slater et al. (2013) show that the southern bifurcated stream is closer to the Sun than the southern portion of the Sgr dwarf tidal stream. Yanny et al. (2009) show that the velocities and metallicities along the bifurcated stream are similar to those in Sgr. In Koposov et al. (2013), there is evidence presented that the two streams may both pass through the progenitor, but due to the proximity of the progenitor to the Galactic bulge, it is difficult to see where exactly the two cross in reference to the progenitor. In Newby et al. (2013), it is suggested the streams may be from two separate progenitors that accreted around the same time, but the evidence to support this is not strong. Hernitschek et al. (2017) give a possible fit to the bifurcated stream. Currently, the origin of this stream is still an open question that our results will help answer. Determining the origin of the "bifurcated" stream is critically important, as it is useful for constraining the Milky Way potential (Law & Majewski 2010; Vera-Ciro & Helmi 2013).

The Virgo Overdensity/Virgo Stellar Stream (Vivas et al. 2001; Duffau et al. 2006; Newberg et al. 2007; Jurić et al. 2008) is a third large halo overdensity in the northern Galactic hemisphere, at distances of 6–20 kpc from the Sun. It is unclear whether this feature is a tidal stream, a "cloud," or a combination of several different pieces. Carlin et al. (2012) fit an orbit to the puffy structure and suggest that this overdensity is the result of a recently disrupted massive (10⁹ M_⊙) dwarf galaxy. Carlin et al. (2012) also find that their orbit includes the Pisces Overdensity. Li et al. (2016b) suggest Virgo could instead be associated with the Hercules-Aquila Cloud and Eridanus-Phoenix overdensities, given that they are on the same polar plane, have similar galactocentric distances (18 kpc), and are separated by 120°. In Bonaca et al. (2012), it is suggested that Virgo is "cloudlike" and may have been the result of a minor merger that passed close to the Galactic center. Vivas et al. (2016) find several different, presumably unrelated, substructures of RR Lyrae stars at distances of 10–20 kpc in the Virgo region. They suggest that there could be additional substructures at much larger distances. The evidence for a more distant Virgo substructure is amplified by Sesar et al. (2017), who find an outer Virgo overdensity at a distance of 80 kpc from the Sun.

The Milky Way stellar halo has traditionally been described by a smooth power-law distribution (e.g., Oort & Plaut 1975; Preston et al. 1991). Since the discovery of significant substructure in the stellar halo (Newberg et al. 2002), researchers have had to choose whether to include or exclude these substructures when fitting the overall spheroid density. The smooth density component of the halo includes smaller or more thoroughly mixed remnants of tidal stripping, as well as any stars that were created in the initial gravitational collapse of the Milky Way galaxy. The algorithms used in this project will fit the smooth component and streams simultaneously. Using simultaneous fitting for the background (as we will refer to the smooth component) and streams, instead of subtracting a background to fit streams, we learn about the shape and density of this smooth component of the stellar halo without requiring a clean sky sample to fit it.

We implemented two different models for the smooth component of the halo, both with two tunable parameters. One model is a Hernquist spheroid model (Hernquist 1990; Xu et al. 2015) with a double exponential disk. The other is a broken power law (BPL) (Akhter et al. 2012) without a disk. In our model fitting, we will primarily use the Hernquist model in order to be consistent with the previous versions of the model described in Cole et al. (2008) and Newby et al. (2013). We first describe the Hernquist/thick disk model, and then the BPL model, which we used to study the effect of an imperfectly modeled smooth component on the derived properties of the halo substructure.

3.1.1. Hernquist Plus Double Exponential Disk

The Hernquist distribution (Hernquist 1990) is described by the equation:

$$\begin{eqnarray}&&{\rho }_{\mathrm{spheroid}}(r)\propto \displaystyle \frac{1}{r{(r+{r}_{0})}^{3}},\end{eqnarray} \tag{ 1 }$$

where $r=\sqrt{{X}^{2}+{Y}^{2}+\tfrac{{Z}^{2}}{{q}^{2}}}$; note that this is not spherical radius, but instead an ellipsoidal radius. In this equation, X, Y and Z are Galactocentric Cartesian coordinates with the Sun at (−8.5, 0, 0) kpc, Y in the direction of the Sun's motion, and Z in the direction of the north Galactic pole. Here, r₀ is a scale radius and q is a flattening parameter. Using this model, there are two tunable parameters: r₀ and q. After inspecting the Hernquist distribution (Hernquist 1990), the Einasto profile (Retana-Montenegro et al. 2012), and a BPL distribution (Akhter et al. 2012), we found that scale radius only has a minor effect on the overall shape of the distribution over the range of our data. In our model, we fix the scale radius at 12.0 kpc, which approximately reproduces the average of all of the distributions inspected.

The double exponential thick disk used in our model is described by the equation:

$$\begin{eqnarray}&&{\rho }_{\mathrm{disk}}({R}_{\mathrm{cyl}},Z)\propto {e}^{{R}_{\mathrm{cyl}}/{l}_{\mathrm{disk}}}{e}^{| Z| /{h}_{\mathrm{disk}}};\end{eqnarray} \tag{ 2 }$$

where ${R}_{\mathrm{cyl}}=\sqrt{{X}^{2}+{Y}^{2}};$ l_disk is the disk scale radius; and h_disk is the disk scale height. In our model, we fix the scale length and height of the disk to be the same as those found in Xu et al. (2015): l_disk = 3500 pc and h_disk = 700 pc. In Figure 2, we show the relative expected fractional contributions of the thick disk, thin disk, and Hernquist background in the Xu et al. (2015) model at different positions in SDSS stripe 19. We do not include the thin disk in our model because our data is far enough from the Galactic plane to avoid the majority of the stars in this component, and our color cuts also help eliminate disk contamination.

Zoom In Zoom Out Reset image size

Because we are using blue MSTO stars that preferentially avoid the thick disk, we cannot use a published normalization for the disk and halo, so this normalization is a fit parameter. We write the overall density of the combined stellar spheroid and disk as:

$$\begin{eqnarray}\begin{array}{rcl}{\rho }_{\mathrm{smooth}}(X,Y,Z) & = & {\epsilon }_{\mathrm{spheroid}}\,\ast \,{\rho }_{\mathrm{spheroid}}(r)\\ & & +\,(1-{\epsilon }_{\mathrm{spheroid}})\,\ast \,{\rho }_{\mathrm{thick}\mathrm{disk}}({R}_{\mathrm{cyl}},Z),\end{array}\end{eqnarray} \tag{ 3 }$$

where _spheroid is a weighting factor that allows us to fit any mix of densities between all disk (_spheroid = 0) and all spheroid (_spheroid = 1). The following equation can be used to convert _spheroid to f_disk(X, Y, Z), the fraction of thick disk stars at a given point:

$$\begin{eqnarray}&&{f}_{\mathrm{disk}}(X,Y,Z)=\displaystyle \frac{(1-{\epsilon }_{\mathrm{spheroid}})\,\ast \,{\rho }_{\mathrm{thick}\mathrm{disk}}({R}_{\mathrm{cyl}},Z)}{{\rho }_{\mathrm{smooth}}(X,Y,Z)}.\end{eqnarray} \tag{ 4 }$$

To find the fraction of stars in the spheroid at a point, f_spheroid(X, Y, Z), use:

$$\begin{eqnarray}&&{f}_{\mathrm{spheroid}}(X,Y,Z)=1-{f}_{\mathrm{disk}}(X,Y,Z).\end{eqnarray} \tag{ 5 }$$

The two parameters that we fit in the smooth Hernquist/double exponential model are q and _spheroid.

3.1.2. BPL Model

The BPL distribution from Akhter et al. (2012) is described by the following equation:

$$\begin{eqnarray}&&{\rho }_{\mathrm{BPL}}(r)\propto {(r/{r}_{0})}^{n},\end{eqnarray} \tag{ 6 }$$

where

$$\begin{eqnarray}&&n=\left\{\begin{array}{l}-2.78\ \mathrm{if}\ r\lt 45\,\mathrm{kpc}\\ -5.0\ \mathrm{if}\ r\geqslant 45\,\mathrm{kpc}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

r₀ = 8.5 kpc is the distance from the Galactic center to the Sun, and $r=\sqrt{{X}^{2}+{Y}^{2}+\tfrac{{Z}^{2}}{{q}^{2}}}$. We use

$$\begin{eqnarray}&&q=\left\{\begin{array}{l}{q}_{0}+(1-{q}_{0})({R}_{\mathrm{cyl}}/{R}_{u})\ \mathrm{if}\ {R}_{\mathrm{cyl}}\lt {R}_{u}\\ 1.0\ \mathrm{if}\ {R}_{\mathrm{cyl}}\gt {R}_{u}\end{array}\right.,\end{eqnarray} \tag{ 8 }$$

where q₀ = 0.5, ${R}_{\mathrm{cyl}}=\sqrt{{X}^{2}+{Y}^{2}}$, and R_u = 20 kpc, as written in Keller et al. (2008).

We do not use this model for fitting the halo because the Hernquist model is more comparable to previous work. Instead, we use this model for testing our algorithm.

In each stripe, which probes a wedge-shaped volume of the Galaxy, the density of a tidal stream or cloud is fit to a cylinder with a uniform density along its length, and Gaussian fall-off in density as a function of radius. The orientation of the stream through the wedge is described by a unit vector $\hat{a}$ with a polar angle from the Z-axis of θ, and an azimuthal angle of ϕ measured from the Galactocentric X-axis and increasing in the direction of the Galactocentric Y-axis. The position of the cylinder's axis, at the point that it passes through the center of the stripe (ν = 0), is given by an angular coordinate along the stripe, μ, and a distance from the Sun, R. The (μ, ν) coordinates are SDSS great circle coordinates York et al. (2000), which are defined separately for each stripe.

At point ${\boldsymbol{p}}$ relative to the cylinder, the stream's density ρ is described by the function:

$$\begin{eqnarray}&&{\rho }_{\mathrm{stream}}({\boldsymbol{p}})\propto {e}^{\tfrac{-{d}^{2}}{2{\sigma }^{2}}},\end{eqnarray} \tag{ 9 }$$

where d is the radial distance from the stream, and σ is the standard deviation describing the stream's width (Cole et al. 2008).

The final parameter in the stream model is the stream weight, ε, which measures the fraction of stars in the wedge that belong to a stream. To get the fraction of stars in a stream or smooth component (spheroid) from the weights, we use the following functions:

$$\begin{eqnarray}\begin{array}{rcl}{f}_{{\mathrm{stream}}_{i}} & = & \displaystyle \frac{{e}^{{\varepsilon }_{i}}}{1+{\displaystyle \sum }_{j=1}^{k}[{e}^{{\varepsilon }_{j}}]}\\ {f}_{\mathrm{spheroid}} & = & \displaystyle \frac{1}{1+{\displaystyle \sum }_{j=1}^{k}[{e}^{{\varepsilon }_{j}}]}.\end{array}\end{eqnarray} \tag{ 10 }$$

Here, i and j denote the stream number for a total of k streams. The weight of a stream is used for model optimization instead of the total star counts because it is a continuous real number that is defined for all real numbers, making it a better parameter than star counts, which are discrete numbers (Cole et al. 2008).

In total, six parameters are required to fit a single stream: the four spatial coordinates, θ, ϕ, μ, and R, which give the stream position and orientation; the width parameter, σ; and the stream weight, ε. The model can handle an arbitrary number of streams in a single wedge.

In our implementation of statistical photometric parallax (Newberg 2013), we use the absolute magnitude distribution for MSTO stars described in Newby et al. (2011), which accounts for the distribution change due to contamination from redder stars at faint magnitudes. In this implementation of statistical photometric parallax, we convolve the model density distributions (in the R direction, which corresponds to our line of sight) with the absolute magnitude distribution. The convolution is given by:

$$\begin{eqnarray}&&\begin{array}{l}{\rho }_{\mathrm{comp}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0}))=\displaystyle \frac{1}{{{ \mathcal R }}^{3}({g}_{0})}{\displaystyle \int }_{-\infty }^{\infty }{{ \mathcal R }}^{3}({g}_{0}-g)\\ \quad \cdot \,{\rho }_{\mathrm{comp}}(l,b,{ \mathcal R }({g}_{0}-g))\cdot N({g}_{0}-g;{g}_{0},u){dg},\end{array}\end{eqnarray} \tag{ 11 }$$

where ${ \mathcal R }({g}_{0})$ is distance to a star (or position in apparent magnitude space) with an apparent magnitude of g₀ and assumed absolute magnitude of ${M}_{{g}_{0}}=4.2$, ρ_comp is a density model, either background or stream, l and b are Galactic l and b, and $N({g}_{0}-g;{g}_{0},u)$ is our MSTO star absolute magnitude distribution. This integral is the combination of the MSTO star absolute magnitude distribution function with our stellar density function over all magnitude space. For an in-depth explanation of this integral, see the original derivation in Cole et al. (2008).

In our model, we use the MSTO absolute magnitude distribution found in Newby et al. (2011), which is defined by two half-Gaussian distributions that are normalized together:

$$\begin{eqnarray}&&N(x;{g}_{0},u)=\displaystyle \frac{1}{\tfrac{({\sigma }_{l}+{\sigma }_{r}({ \mathcal R }({g}_{0})))}{2}\sqrt{2\pi }}{e}^{\tfrac{-{x}^{2}}{2{u}^{2}}},\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&u=\left\{\begin{array}{l}{\sigma }_{l}=.36\ \mathrm{if}\ x\lt 0\\ {\sigma }_{r}({ \mathcal R }({g}_{0}))=\displaystyle \frac{\alpha }{1+{e}^{-({ \mathcal R }({g}_{0})-\beta )}}+\gamma \ \mathrm{if}\ x\geqslant 0\end{array}\right.,\end{eqnarray} \tag{ 13 }$$

α = 0.52, β = 12.0, and γ = 0.76. This absolute magnitude distribution accounts for an increase in the standard deviation on the faint side of the absolute magnitude distribution as a function of magnitude in the SDSS. The change in distribution is due to the change in the type of stars found in our narrow MSTO color selection criteria as increasing color errors toward the magnitude limit of the survey. Toward the survey magnitude limit, color errors increase and cause redder, lower main sequence stars to leak into our color selection as actual MSTO stars leak out. The underlying absolute magnitude distributions of all turnoff stars in the halo is assumed to be constant as found by Newby et al. (2011).

To correct for the drop-off in detection efficiency near the magnitude limit of the SDSS, a sigmoid curve was fit in Cole (2009) to data from Newberg et al. (2002):

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{sigmoid}}({g}_{0})=\displaystyle \frac{{s}_{0}}{{e}^{{s}_{1}({g}_{0}-{s}_{2})}+1},\end{eqnarray} \tag{ 14 }$$

where (s₀, s₁, s₂) = (0.9402, 1.6171, 23.5877). This curve describes the survey completeness at a given magnitude. Using this curve, we can account for the decreasing completeness of the survey with increasing magnitude.

As the color errors increase, fainter main sequence stars leak into the MSTO color selection box, and MSTO stars leak out of the MSTO color selection box. The detection efficiency is adjusted using what we call the selection efficiency to account for the number of stars expected given this effect. The selection efficiency for SDSS stars in the MSTO selection bin is given by a seventh-order polynomial fit by Newby et al. (2011). We reproduce it here because there was an inadvertent truncation of significant digits in the original paper that affects the result of the fit:

$$\begin{eqnarray}&&{{ \mathcal E }}_{\mathrm{selction}}({ \mathcal R }({g}_{0}))=\displaystyle \frac{n({ \mathcal R }({g}_{0}))}{{n}_{0}}=\displaystyle \sum _{i=0}^{7}({a}_{{yi}}+{a}_{{ri}}){ \mathcal R }{({g}_{0})}^{i},\end{eqnarray} \tag{ 15 }$$

where a_y = (1.05628761, −3.14555041 × 10⁻², 2.05499665 × 10⁻⁴, 2.53747387 × 10⁻⁶, −2.67000303 × 10⁻⁸, 0, 0, 0) and represents the remaining MSTO stars in the color selection bin, and a_r = (1.60879353 × 10⁻², −1.97164570 × 10⁻², −4.31844102 × 10⁻⁴, 6.60960070 × 10⁻³, 1.26368065 × 10⁻⁵, −1.91560491 × 10⁻⁷, 1.47140445 × 10⁻⁹, −4.53857248 × 10⁻¹²) and represents the red stars that leak into the color selection bin to contaminate the sample. The percentage of stars left in the selection bin, assuming 100% detection efficiency, is illustrated in Figure 3. The two detection efficiency functions are multiplied together to get the total detection efficiency for stars that fit the MSTO selection criteria.

Zoom In Zoom Out Reset image size

Putting together all of the individual model components, we create an MLE that tells us how well our model, with parameters ${\boldsymbol{Q}}$, fits a set of data (N stars with measured l_i, b_i, g_i):

$$\begin{eqnarray}&&{ \mathcal L }({\boldsymbol{Q}})=\displaystyle \prod _{i=1}^{N}\mathrm{PDF}({l}_{i},{b}_{i},R({g}_{i})| {\boldsymbol{Q}}),\end{eqnarray} \tag{ 16 }$$

where i is the index of each star, and l and b are Galactic coordinates (Cole et al. 2008). Substituting in our density functions, correcting for imperfect knowledge of the absolute magnitudes of the stars, and taking into account the detection efficiency, our PDF is given by:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{PDF}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})\\ =\,\displaystyle \frac{1}{1+{\displaystyle \sum }_{i=1}^{n}{e}^{{\epsilon }_{i}}}\displaystyle \frac{{ \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{\mathrm{background}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})}{\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{\mathrm{background}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){dV}}\\ +\,\displaystyle \sum _{i\,=\,1}^{n}\displaystyle \frac{{e}^{{\epsilon }_{i}}}{1+{\displaystyle \sum }_{j\,=\,1}^{n}{e}^{{\epsilon }_{j}}}\displaystyle \frac{{ \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{{\mathrm{stream}}_{i}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})}{\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{{\mathrm{stream}}_{i}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){dV}},\end{array}\end{eqnarray} \tag{ 17 }$$

where n is the number of streams in the model and ${ \mathcal E }({ \mathcal R }({g}_{0}))={{ \mathcal E }}_{\mathrm{sigmoid}}({ \mathcal R }({g}_{0}))\cdot {{ \mathcal E }}_{\mathrm{selection}}({ \mathcal R }({g}_{0}))$. To avoid problems with numerical underflow due to small probabilities, we instead maximize the log likelihood:

$$\begin{eqnarray}&&\displaystyle \frac{1}{N}\mathrm{ln}({ \mathcal L }({\boldsymbol{Q}}))=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}\mathrm{ln}(\mathrm{PDF}({l}_{i},{b}_{i},{g}_{i}| {\boldsymbol{Q}})).\end{eqnarray} \tag{ 18 }$$

In the SDSS, there are areas of the sky with either missing or unusable data, such as the regions hidden behind a saturated bright star. There are also other sections of the sky with known structures that are not included in our model, such as globular clusters. We remove regions with these substructures or artifacts from our data and then integrate over the remaining usable survey volume.

The technical details of this process are as follows. We change our PDF to remove the integral of the density over the volume we are cutting from the denominator of that density's respective fraction:

$$\begin{eqnarray}&&\begin{array}{l}\mathrm{PDF}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})\\ =\,\displaystyle \frac{1}{1+{\displaystyle \sum }_{i\,=\,1}^{n}{e}^{{\epsilon }_{i}}}\displaystyle \frac{{ \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{\mathrm{background}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})}{\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{\mathrm{background}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){dV}-\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{\mathrm{background}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){{dV}}_{\mathrm{cut}}}\\ +\,\displaystyle \sum _{i\,=\,1}^{n}\displaystyle \frac{{e}^{{\epsilon }_{i}}}{1+{\displaystyle \sum }_{j\,=\,1}^{n}{e}^{{\epsilon }_{j}}}\displaystyle \frac{{ \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{{\mathrm{stream}}_{i}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}})}{\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{{\mathrm{stream}}_{i}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){dV}-\displaystyle \int { \mathcal E }({ \mathcal R }({g}_{0})){\rho }_{{\mathrm{stream}}_{i}}^{\mathrm{con}}(l,b,{ \mathcal R }({g}_{0})| {\boldsymbol{Q}}){{dV}}_{\mathrm{cut}}},\end{array}\end{eqnarray} \tag{ 19 }$$

where dV_cut represents integrating over the volume that was cut out of the data and the rest is the same as Equation (17).

In our data, we will either remove the stars from the globular cluster from the data set, or we will not have data in this region due to deblending limitations in the SDSS survey.

There are several degeneracies and covariances we must account for in our parametric model. The parameters that are explicitly covariant are the relative weights between the density of the streams and background. There are many other covariances that can be read from the error matrix, in Table 1. Matrix elements with a high value have high covariance in the corresponding parameters (Ivezić et al. 2014). When comparing two parameters, the sign of the covariance tells you the orientation of the error ellipse. A positive covariance indicates they are positively correlated and a negative covariance indicates they are inversely correlated. An example of parameters with a high covariance is the angular position along the stripe, μ, and the distance from the Sun, R, for the streams.

The degeneracy in stream orientation and stream selection are inherent to an optimization problem with a model of this type. The stream orientation degeneracy comes about due to the cyclic nature of angles and can be mitigated by constraining the optimization to a single hemisphere. The stream selection degeneracy arises from the freedom to fit each stream to any location in the wedge. Therefore, one optimization may find stream 1 to best fit Sagittarius, while another optimization on the same data may find stream 2 to best fit Sagittarius with the same parameters. In this instance, we do not try to eliminate this degeneracy because it does not change the result in any substantive way, and attempting to force a particular stream order would require us to artificially constrain streams to mutually exclusive regions and could bias our optimization.

To test that both the model and fitting algorithm work as intended, the algorithm was tested on a set of data drawn from a density model with known parameters. We wrote a completely separate program to generate the test data, so that the entire modeling and fitting procedure would be tested. Generating the test data is a multistep process in which we simulate our background stellar distribution and streams with the correct number of stars, determine the star's absolute and apparent magnitudes using our magnitude distribution, and then apply the effects of observational bias.

The first step in generating our test data is selecting the set of parameters we wish to simulate. To do this, we use the values of the parameters from our fit to SDSS stripe 19 from Table 2. This is the best set of parameters to check the accuracy of the derived stripe 19 density model.

After we choose the parameters, we input them to our wedge simulator, which begins by calculating the number of stars in each component of our data. The program does this using Equation (10) with the specified weights to get the fraction of stars in each component. The fraction of stars is then multiplied by the number of MSTO stars we want to simulate (in this case, 84,046, which is the number of turnoff stars observed in stripe 19). Each component is then simulated with the required number of stars.

For each component of the density, the test data generator creates one star at a time until the number required by the weights is reached. The smooth component is generated using rejection sampling. We uniformly sample Equation (3) over the volume of the wedge, and we compare the density from Equation (3) to the maximum density possible in the wedge from Equation (3). Then, using a random number, we determine whether we will reject or keep the star, based on the fraction of the calculated density to the maximum density.

Stream stars are generated using active generation to create a stellar position in three dimensions. One number drawn from a uniform distribution determines the position along the cylinder's z-axis, and two numbers drawn from a normal distribution determine the position of the star in the plane perpendicular to the cylinder's z-axis. The length of the cylindrical distribution that is generated is significantly longer than the portion of the cylinder axis that is within the wedge of data in stripe 19. The cylinder's length is long enough that all of the volume in the wedge that is within three σ of the cylinder's z-axis will be populated. The stream star is then rotated and translated into Galactic X, Y, Z based on the stream's orientation and center point.

As the program generates each star, it applies the observational biases that are corrected for in our maximum likelihood model. First, it determines the star's absolute magnitude based on the absolute magnitude distribution in Equation (13) and then uses that to determine its apparent magnitude. At this point, it checks to make sure the star lies within the wedge. Finally, it samples the combined magnitude limit from Equation (14) and the color selection efficiency from Equation (15) to determine if the star would have been observed.

We generated two different sets of test data: one with a Hernquist background that exactly matches the model to which the data is fit, and one with a BPL background. The second simulation tests the sensitivity of our measurements of the stream parameters to imperfect knowledge of the smooth component of the Milky Way spheroid. For the BPL background, we used the parameters fit in Akhter et al. (2012): an inner halo exponent of −2.78, outer halo exponent of −5.0, and break radius of 45 kpc.

MilkyWay@home was used to find the optimum values of the parameters, fit separately for the two sets of generated test data. We ran four optimizations for each of the simulated data wedges, each with minimal constraints. Each set of four optimizations independently converged to the same results; these results are listed in Table 4.

Next, we determined the uncertainty in each parameter, using the method outlined in Section 5. Because our errors are given as one standard deviation, we can expect 68% (13–14 parameters) of the results to be within uncertainty, 95% (19 parameters) to be within two times the uncertainty, and the rest to be within three times the uncertainty, most of the time. Our results in Table 4 show that, for the Hernquist plus disk background simulation, we had 10 parameters within the uncertainties, 16 within two times the uncertainties, and 18 were within three times the uncertainties. While it is unusual that two of the parameters are more than three times the uncertainty from their simulated value, we do not believe this is enough to cause alarm.

In our results for the simulation with the BPL model, we find that five parameters are within the uncertainties, eight are within two times the uncertainties, and nine are within three times the uncertainties. The higher frequency of best-fit parameters with large offsets from the simulated values shows that, if our model is not a good match to the actual stellar density model in the halo, our errors, as calculated using the method in Section 5, could be underestimated by a factor of three or more. The fraction of stars in each component is particularly poorly measured compared to the uncertainties; the optimizations for ε differ from the simulated values by 6–11σ. These errors look worse than they really are because they are measured with respect to a different smooth component. We include a column in Table 4 that shows the number of stars in the stripe associated with each component. The density of stars along a stream could be off by as much as 60%.

Nevertheless, the properties of the detected streams and the general shape of the stellar spheroid (as measured by q) are generally similar to the simulation. All of the halo substructure is immediately recognizable. The flattening of the spheroid is surprisingly correct. The widths of the streams are approximately correct, the stream centers and distances are pretty close, and even the angles at which the steams pass through the wedge are easily matched with the correct stream.

In addition to checking the accuracy of the uncertainties, we checked the accuracy of the optimizer in finding the highest point in the likelihood surface. To do this, we ran one-dimensional parameter sweeps for each of the 20 parameters in our model. The sweeps seen in Figure 6 show the likelihood surface for the test data with the Hernquist background. This figure shows that each model parameter peaks either at, or close to, the expected parameter value. Any small difference between the peak in the likelihood surface and the simulated values can be accounted for by the uncertainty introduced by the finite number of stars available in our data. For the BPL parameter sweeps in Figure 7, there are several parameters that have extremely flat likelihood surfaces. In these cases, the optimizer still seems to find an answer close to the peak of the likelihood.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to visually compare the results of our simulation to the actual SDSS stripe 19 data (Figure 4), we present the simulated data with a Hernquist background in Figure 8, and the simulated data with the BPL background in Figure 9. The simulated data figures (Figures 8 and 9) show three different views of the data. The first row shows the components generated by the test data simulator, and the combination of those components. The second row shows a probabilistic separation of the components by our model, given the simulation parameters. Finally, the last row shows a resimulation of the data using the fit model parameters. By resimulating the stripe, we can visualize the model that the optimizer found as the best fit to the data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In each of these plots, the separated streams have the appearance of streams, and the smooth component appropriately contains no substructure. It is especially interesting that MilkyWay@home is able to determine, on its own, a set of parameters that resulted in streams that look, by eye, the same as those with the "correct" parameters for the wedge. More evidence of this can be seen in Figures 10 and 11, which show residual densities of the simulated data and resimulated fit data from Figures 8 and 9. The BPL fits show that we will not easily know whether or not the real data is a poor fit to a Hernquist background. However, we will still be able to fit most stream parameters (eight of the 18 stream parameters were within 2σ of the simulated values), even though a few might not be correct within the calculated errors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In summary, simulations show that in the case that our background and stream models are close to the correct answer, the resultant parameter fits can be trusted. If our background or stream models are far from the correct distribution of stars, the uncertainties are underestimated, but the stream properties are generally reasonable.