TOUCHING THE VOID: A STRIKING DROP IN STELLAR HALO DENSITY BEYOND 50 kpc

In our model universe, the balance between expansion and collapse stipulates that the size and the mass of a galaxy are set by its formation epoch (see, e.g., Press & Schechter 1974). Once most of the galactic contents are in place, subsequent matter infall adds little to the final mass budget (see, e.g., Zemp 2013). The total mass is dominated by dark matter; even though gas and stars might extend as far, their densities drop faster with radius and therefore contribute little to the integral over the virial volume. However, despite amounting to only 1% of the total galaxy luminosity or <0.01% of the total mass, the stellar halo allows us to gauge the details of the mass distribution beyond the edge of the disk.

The stars in the halo are more than mere tracers of the potential. The dark matter radial density profiles are universal (and hence featureless beyond the scale radius), or at least they appear to be so for a considerable range of distances explored in numerical simulations (e.g., Navarro et al. 1997). However, due to the plummeting star-formation efficiency in low-mass sub-halos (e.g., Bullock et al. 2000; Somerville 2002), the stellar halo formation is a much more stochastic process. The lumpier accretion, combined with extremely long mixing times (>1 Gyr) can lead to a greater variety of stellar halo radial density profiles (see, e.g., Libeskind et al. 2011). Therefore, there is hope that by studying the phase-space and chemical properties of halo stars today, we can uncover the fossil record of the Milky Way's accretion history.

In order to quantify the stellar halo distribution, we often fit model profiles, such as power-laws and Einasto profiles (Einasto & Haud 1989), to the stellar number counts. This approach has been widely used in the literature, and although these models may not represent a truly physical representation of the stellar halo, they provide a useful framework that can be compared with predictions from numerical simulations. Early work limited to Galactocentric radii r ∼ 20–30 kpc found that the Milky Way stellar halo follows an oblate, single power-law distribution with minor-to-major axis ratio q ∼ 0.5–0.8, and power-law index α ∼ 2–4 (e.g., Preston et al. 1991; Robin et al. 2000; Yanny et al. 2000; Newberg & Yanny 2006; Jurić et al. 2008). More recent work, probing to greater distances in the halo, found evidence for a "break" in the stellar density profile at r ∼ 20–30 kpc.⁵ These studies find a power-law slope of α ∼ 2–3 can describe the stellar halo within r ∼ 20–30 kpc, but a steeper slope with α ∼ 3.8–5 is required at larger distances (Bell et al. 2008; Watkins et al. 2009; Sesar et al. 2011; Deason et al. 2011a).

Deason et al. (2013a) argued that this broken profile could be caused by the build-up of stars at their apocenters, either from the accretion of one massive dwarf, or from several dwarfs accreted at a similar epoch. On the other hand, Beers et al. (2012) claim that the change in power-law slope near the break radius is caused by a transition from an "inner" to an "outer" stellar halo population. Several groups have found evidence for correlations between metallicity and kinematics of halo stars, which perhaps suggest two distinct populations (e.g., Carollo et al. 2007, 2010; Nissen & Schuster 2010; Deason et al. 2011b; Hattori et al. 2013; Kafle et al. 2013). However, at present it is not obvious whether these signatures can be produced purely from the accretion of dwarf galaxies, or if some of these findings are biased by distance uncertainties and/or contamination (e.g., Schönrich et al. 2011, 2014; Fermani & Schönrich 2013).

It is clear that from a relatively "simple" measure of star counts, we can learn a great deal about the formation mechanism and/or past accretion history of the stellar halo. This bodes well for studies of stellar halos beyond the local group, where we are already able to measure the surface brightness profiles of these incredibly diffuse halos out to projected radii of R ∼ 50–70 kpc (e.g., Radburn-Smith et al. 2011; Monachesi et al. 2013; Greggio et al. 2014; van Dokkum et al. 2014). The surface brightness profile of our nearest neighbor, M31, has now been mapped out to an impressive R ∼ 200 kpc (e.g., Gilbert et al. 2012; Ibata et al. 2014). In contrast to our own Galaxy, these studies find no evidence for a break in the stellar density profile, and the star counts can be well-described by a single power-law with slope α ∼ 3–3.5. The differences between the stellar halo density profiles can give us an important insight into the contrasting accretion histories of the Milky Way and M31 (see Deason et al. 2013a).

Superimposed on the "field" (or phase-mixed) stellar halo distribution is a wealth of un-relaxed substructure in the form of streams, clouds and other overdensities (e.g., Ibata et al. 1995; Newberg et al. 2002; Belokurov et al. 2006; Belokurov et al. 2007; Jurić et al. 2008). Most striking is the vast stream of tidal debris associated with the disrupting Sagittarius dwarf (Belokurov et al. 2006), where we are privy to a front-seat view of accretion in action. The presence of these recent accretion relics can significantly bias stellar halo number counts. Several studies have attempted to excise these known substructures, and only model the relatively phase-mixed halo component. However, it is important to understand the affect that these structures have on density profile measurements, especially when making comparisons with numerical simulations or stellar halos of external galaxies. For the latter, it is generally unfeasible to isolate the sort of substructures that we are able to identify in our own Galaxy.

Halo stars provide one of the best tracers of the Milky Way mass at large radii. Deason et al. (2012) compiled a sample of stellar halo stars with measured line-of-sight (LOS) velocities out to r ∼ 150 kpc, and found a dramatic drop in LOS velocity dispersion beyond 50 kpc (see also Battaglia et al. 2005). If this drop reflects a fall in the circular velocity of the halo, then the Milky Way mass is likely below ∼10¹² M_☉. This result agrees with several other stellar dynamical studies, which favor relatively low halo masses (e.g., Xue et al. 2008; Bovy et al. 2012). By contrast, other methods for estimating the Milky Way mass give larger values. For example, Sohn et al. (2013) recently used multi-epoch Hubble Space Telescope (HST) images to measure the proper motion of the distant Milky Way satellite galaxy Leo I. Comparison of the large observed velocity to numerical simulations implies that the Milky Way mass is likely well above 10¹² M_☉ (Boylan-Kolchin et al. 2013). More generally, attempts to measure the total mass of the Milky Way using satellite galaxies (e.g., Wilkinson & Evans 1999; Watkins et al. 2010), the Magellanic Clouds (e.g., Kallivayalil et al. 2013), the local escape speed (e.g., Smith et al. 2007), the timing argument (e.g., Li & White 2008; van der Marel et al. 2012; Gonzalez et al. 2013), and the methods already mentioned, have been distressingly inconclusive with total masses in the range 0.5–3 × 10¹² M_☉.

Halo stars have tremendous potential for constraining the Milky Way mass, since LOS velocities have been measured for many of them, but to make progress the mass–anisotropy–density degeneracy must be addressed. Our mass measures based on halo star kinematics are limited by the uncertainty in the tracer density profile and velocity anisotropy. These systematic uncertainties are significant, and mass-measures can vary by up to factors of ∼5 because of unknown tracer properties. Fortunately, the upcoming Gaia mission and deep, multi-epoch HST proper motion measurements (Deason et al. 2013b; van der Marel et al. 2013), will provide the missing transverse velocity information needed to measure the velocity ellipsoid of distant halo stars. However, we still have very little knowledge of the tracer density profile beyond 50 kpc. Thus, somewhat ironically, the "simple" task of counting stars will likely be the main bottleneck for dynamical mass measures of the Milky Way in the near future.

In this study, we use A-type halo stars selected from Sloan Digital Sky Survey (SDSS) data release 9 (DR9) photometry to measure the stellar halo density slope beyond ∼50 kpc. These A-type stars comprise of blue horizontal branch (BHB) and blue straggler (BS) populations. The former stellar population, constitute our prime halo tracers and can probe out to ∼100 kpc in the magnitude range used in this work (18.5 < g < 20.5). Our method models both BHB and BS populations simultaneously using photometric data alone, and includes the contribution from contaminants, such as QSOs. The combination of the large SDSS sky coverage (∼14, 000 deg²) and the accurate distance estimates provided by the BHB stars, allows us, for the first time, to constrain the outer density profile slope of the Milky Way stellar halo.

The paper is arranged as follows. In Section 2.1 we describe the SDSS DR9 photometric data and our selection criteria for A-type stars. The remainder of Section 2 describes our A-type star models and the absolute magnitude–color relations for the two populations. In Section 3 we address the contribution of contaminants and the affects of photometric uncertainties on our modeling procedure. In Section 4, we describe our likelihood-based method to determine the density profile of the stellar halo and in Section 5 we present our results. Finally, we discuss the implications for the accretion history and the mass of the Milky Way in Section 6, and summarize our main conclusions in Section 7.

2.1. DR9 Imaging

The SDSS (York et al. 2000) is an imaging and spectroscopic survey covering over one quarter of the sky. Images are obtained simultaneously in five broad optical bands (ugriz; Fukugita et al. 1996) using a CCD camera (Gunn et al. 1998) on a 2.5 m telescope (Gunn et al. 2006) at Apache Point Observatory, New Mexico. Photometric and astrometric properties are derived through data processing pipelines developed throughout the course of the survey (Lupton et al. 2001; Smith et al. 2002; Stoughton et al. 2002; Pier et al. 2003; Ivezić et al. 2004; Tucker et al. 2006). The SDSS DR9 (Ahn et al. 2012) provides the same sky coverage as its predecessor DR8 (∼14, 000 deg²), and contains all of the imaging data taken by the SDSS imaging camera.⁶ We select high latitude (|b| > 30 deg) objects classified as stars by SDSS with clean r-band photometry. The magnitudes and colors we use in the following sections have been corrected for extinction following the prescription of Schlegel et al. (1998).

In this study, we use BHB and BS A-type stars to map the density profile of the distant Milky Way stellar halo. These high latitude A-type stars are limited to a tight sequence in u − g, g − r space (see text below and Figure 1) and suffer from relatively little contamination. This, in combination with the well-determined photometric parallaxes of BHB stars (see Section 2.2.1), make these ideal tracers of the Milky Way halo. Deason et al. (2011a; hereafter, DBE11) used photometrically selected A-type stars from SDSS DR8, with 16.0 < g < 18.5, to determine the stellar halo density profile within r ∼ 40–50 kpc. Here, we extend this analysis to fainter magnitudes, 18.5 < g < 20.5, in order to measure the stellar density profile beyond 50 kpc. In the magnitude range under consideration, BS/BHB stars have approximate heliocentric distances: 10 ≲ D_BS/kpc ≲ 75, 40 ≲ D_BHB/kpc ≲ 100.

Zoom In Zoom Out Reset image size

In Figure 1 we indicate our A-type star selection in u − g, g − r (or ugr for short) space, 0.7 < u − g < 1.6, −0.25 < g − r < −0.1, as a function of g-band magnitude. The range in u − g is slightly expanded relative to DBE11 (cf. 0.9 < u − g < 1.4) to account for the increased photometric errors at fainter magnitudes (see Figure 2). In addition, we impose a bluer g − r selection (g − r < −0.1, DBE11 used g − r < 0) to minimize contributions from redder BS stars (see e.g., Brown et al. 2010) and QSOs. At fainter magnitudes the A-type "claw" sequence at u − g ∼ 1 widens and the broad QSO distribution begins to impose on our A-type star selection box (see Section 3.1.1).

Zoom In Zoom Out Reset image size

2.2. A-type Star Models

The modeling procedure used by DBE11 assumed that only A-type stars (BHB or BS) were present in their sample, and they did not suffer from significant photometric errors. Our extension to fainter magnitudes means that we need to consider contamination from other populations, and the affect of photometric scattering; these model extensions are discussed in Sections 3.1 and 3.2.

Our A-type star models describe the distribution of BHB and BS stars in ugr, m_g space. The three-dimensional model for each stellar population depends on (1) the intrinsic ugr color distribution, (2) the stellar density profile and (3) the absolute magnitude-(g − r) color relation.

Owing to the higher surface gravity of BS stars, BHB and BS populations can be distinguished spectroscopically from their Balmer line profiles (e.g., Kinman et al. 1994; Clewley et al. 2002; Sirko et al. 2004; Xue et al. 2008; Deason et al. 2012). However, with photometry alone BHB and BS stars form distinct, but overlapping sequences in u − g, g − r color–color space. DBE11 used bright A-type stars with available SDSS spectra to pinpoint the loci of BHB and BS populations in u − g, g − r color–color space. These "ridgelines," i.e., the approximate centers of their distributions in u − g as a function of g − r, are defined by third order polynomials:

$$\begin{eqnarray} (u-g)^{0}_{\rm BHB}&=& 1.167-0.775(g-r)-1.934(g-r)^2 \nonumber \\ &&+9.936(g-r)^3, \nonumber \\ (u-g)^0_{\rm BS}&=& 1.078-0.489(g-r)+0.556(g-r)^2 \nonumber \\ &&+13.444(g-r)^3, \end{eqnarray} \tag{ 1 }$$

for −0.25 < g − r < 0.0. (see left-hand panel of Figure 1, and DBE11, Figure 2)

The intrinsic spread of the two populations about their ridgelines are σ_{BHB, 0}(u − g) = 0.04 and σ_{BS, 0}(u − g) = 0.045. Gaussian distributions about these ridgelines are adopted to assign a membership probability based on u − g color:

$$\begin{eqnarray} && p(u-g|{\rm BHB}, g-r)\propto \mathrm{exp}\left(-\frac{\left[(u-g)-(u-g)_{\rm BHB}^0\right]^2}{2\sigma _{\rm BHB}^2}\right), \nonumber \\ && p(u-g|{\rm BS}, g-r)\propto \mathrm{exp}\left(-\frac{\left[(u-g)-(u-g)_{\rm BS}^0\right]^2}{2\sigma _{\rm BS}^2}\right). \end{eqnarray} \tag{ 2 }$$

Here, the ridgelines (see Equation (1)), (u − g)⁰, depend on g − r color.

The g − r distribution of A-type stars depends on the stellar density profile, and the intrinsic g − r distribution of each population, i.e., the relative volume densities of objects with different g − r. Thus, without a priori knowledge of the stellar density profile, we cannot disentangle the intrinsic g − r distribution. It is important that the apparent g − r distribution for BS stars (and to a lesser extent BHB stars) at fixed g-band magnitude is different from the intrinsic one, and is a (strong) function of the density profile.

DBE11 estimated the relative fractions of BHB and BS stars in g − r bins using the u − g model distributions defined above. However, in the fainter magnitude range under consideration here, the u − g random errors are too large to be able to adequately distinguish between BHB and BS stars (see Figure 2 where σ_{(u − g)} ∼ 0.05–0.1 between 18.5 < g < 20.5). Instead, we use the stellar density profile derived by DBE11 for bright A-type stars (16 < g < 18.5), and the approximate numbers of BHB and BS stars in g − r bins (see Table 1 in DBE11) to disentangle the intrinsic g − r distributions of these populations. Note that this exercise assumes that the intrinsic distributions remain constant with magnitude.

In Figure 3 we show the resulting intrinsic g − r distributions for (bright) BHB and BS stars. The BS stars have a steep g − r distribution which rises sharply toward redder g − r color, while the BHB stars are roughly constant with g − r. The dotted lines indicate polynomial fits to the g − r distributions: p(g − r|BHB)∝0.70 − 44.7(g − r) − 144.7(g − r)², p(g − r|BS)∝26.6 + 166.8(g − r) + 265.7(g − r)². Note that while the intrinsic g − r distributions of BHB/BS stars are an important ingredient of the modeling procedure, our results are not significantly affected by our adopted parameterization. For example, our main conclusions are unchanged if we instead adopt a flat g − r distribution for both populations.

Zoom In Zoom Out Reset image size

Our A-type star membership probabilities can now be assigned based on u − g, g − r colors:

$$\begin{eqnarray} &&p({ugr}|{\rm BHB})\propto p(g-r | \mathrm{BHB}) \, p(u-g | \mathrm{BHB}, g-r), \nonumber \\ &&p({ugr}|{\rm BS})\propto p(g-r | \mathrm{BS}) \, p(u-g | \mathrm{BS}, g-r). \end{eqnarray} \tag{ 3 }$$

2.2.1. Absolute Magnitude Calibration

BHB stars are intrinsically brighter than BS stars (by ∼2 mag), and their absolute magnitude varies little as a function of temperature or metallicity. By comparison, BS stars are intrinsically fainter and span a much wider range in absolute magnitude. We adopt the DBE11 absolute magnitude–color relations for these populations. DBE11 used star clusters with SDSS photometry published by An et al. (2008) to calibrate the BHB absolute magnitudes, and stars selected in the Sagittarius stream with Stripe 82 photometry were used to calibrate the BS star absolute magnitudes (see Watkins et al. 2009). The derived relations are repeated here for completeness:

$$\begin{eqnarray} M_{g(\rm BHB)}&=& 0.434-0.169(g-r)+2.319(g-r)^2 \nonumber \\ &&+20.449(g-r)^3+94.517(g-r)^4, \nonumber \\ M_{g(\rm BS)}&=& 3.108+5.495(g-r), \end{eqnarray} \tag{ 4 }$$

where, $\sigma _{M_g(\rm BS)} \sim 0.5$.

3.1. Contamination

Here, we consider the affect of contaminants on our A-type star selection box. In the left-hand panel of Figure 4 we show the u − g distribution for objects with −0.25 < g − r < −0.1 and 18.5 < g < 20.5. The purple dashed lines indicate our u − g selection boundary. QSOs and white dwarfs (WDs) have bluer u − g colors than A-type stars, and these populations can clearly be seen for u − g < 0.6. Due to the relatively large photometric errors in the magnitude range under consideration, we must take into account these populations in our modeling procedure.

Zoom In Zoom Out Reset image size

The blue dotted line shows the u − g distribution when a cut in gri color space is applied to exclude QSOs (see Deason et al. 2012, Figure 2). This cut is able to remove a significant amount of QSOs, but, owing to the large photometric errors at fainter magnitudes, it cannot remove all of them. Furthermore, this cut would also remove some A-type stars, especially those at fainter magnitudes with larger photometric errors. Thus, in order to avoid any biases, we do not apply this cut, but instead include a model for the QSO population in our analysis. The difference between the u − g distributions with and without the gri QSO cut illustrates that the QSOs constitute our main contaminant. In the following section we construct a model for the QSO population.

3.1.1. QSO Models

To model the QSO population we make use of the XDQSO code⁷ which is designed to calculate photometric QSO probabilities (see Bovy et al. 2011a). The XDQSO algorithm was developed by Bovy et al. (2011a) for efficient flux-based QSO target selection for SDSS data. Models of quasars in flux space were built by applying the extreme-deconvolution method (Bovy et al. 2011b) to spectroscopically confirmed QSOs in order to estimate the underlying density. This density is convolved with the flux uncertainties when evaluating the probability that an object is a QSO. We use these intrinsic photometric QSO models to construct the unconvolved (i.e., not convolved with photometric errors) QSO probability density function (PDF) in ugr, m_g space.

The XDQSO code is applied to a uniform ugriz distribution in the appropriate color and magnitude range. The QSO PDF is generated from the output likelihood distributions for low, medium and high redshift quasars:

$$\begin{equation} P\mathrm{(QSO)}=\mathcal {L}_{z_{\rm low}}N_{z_{\rm low}}+\mathcal {L}_{z_{\rm mid}}N_{z_{\rm mid}}+\mathcal {L}_{z_{\rm high}} N_{z_{\rm high}}. \end{equation} \tag{ 5 }$$

Here, $\mathcal {L}_{z_{\rm low}}$, $\mathcal {L}_{z_{\rm mid}}$, and $\mathcal {L}_{z_{\rm high}}$ are the relative flux likelihoods for low, medium, and high redshift QSOs, and $N_{z_{\rm low}}$, $N_{z_{\rm mid}}$, and $N_{z_{\rm high}}$ are the number counts at a given i-magnitude. The XDQSO algorithm computes the likelihoods and number counts for each QSO class given ugriz fluxes. An acception–rejection algorithm is used to draw u − g, g − r, and m_g values from this PDF, which gives the unconvolved QSO PDF in ugr, m_g space:

$$\begin{equation} \nu _Q=p(ugr, m_g | \rm QSO). \end{equation} \tag{ 6 }$$

The resulting u − g, g − r color distribution of the unconvolved QSO PDF is shown in the left-hand panel of Figure 5. The purple box indicates our A-type star selection box. The inclusion of photometric errors scatters a significant number of QSOs into the box. The g-band magnitude distribution is shown in the right-hand panel. The QSO luminosity function increases steeply with magnitude, thus the QSO contribution becomes more significant at fainter magnitudes.

Zoom In Zoom Out Reset image size

3.1.2. White Dwarf Models

We now consider the influence of the WD population in our A-type star selection box. In the right-hand panel of Figure 4 we show the u − g distribution for low- and high-probability QSOs determined from the XDQSO algorithm. The low-probability QSOs at bluer u − g are WD stars. The blue dashed line shows the predicted u − g distribution of our WD model (outlined in Appendix A), where suitable photometric errors (see Figure 2) have been included. Our WD model predicts a very small fraction of WD stars in our A-type star selection box (<1%). Thus, we can safely ignore the contribution of WD stars for the remainder of our analysis.

3.2. Photometric Scattering

Our A-type star and QSO models describe the intrinsic, unconvolved populations. In the following section, we apply our likelihood method to the convolved models, which takes into account the ugr photometric uncertainties and their dependence on g-band magnitude. To illustrate the importance of taking photometric errors into account, we show in Figure 6 an estimate of the completeness of our SDSS DR9 catalog. To estimate the completeness of our A-type sample we make use of the ∼2 mag deeper stacked Stripe 82 photometry (Annis et al. 2011). Stars in the color range (0.7 < u − g < 1.6, −0.25 < g − r < −0.1) are selected from Stripe 82 and cross-matched with our DR9 catalog. The black line indicates the completeness fraction based on Stripe 82 stars which are detected in the DR9 catalog. The red line shows the completeness fraction when we also impose that the same stars fall within our ugr selection box. Purely based on detection (i.e., the existence of an object with particular coordinates in both catalogs), our catalog is close to 95% complete over the magnitude range under consideration (18.5 < g < 20.5). However, when we impose that the stars in Stripe 82 are also detected in the same ugr range in the DR9 data, our recovered fraction is reduced, even at brighter magnitudes. This is due to photometric scattering; stars can be scattered into and out of our selection box and it is important to take this into account, especially at fainter magnitudes.

Zoom In Zoom Out Reset image size

The affect of photometric scattering can depend strongly on the shape of the PDF. For example, the intrinsic distribution of BS stars depends strongly on g − r color, whereas BHB stars have a much weaker dependence on g − r (see Figure 3). The steep gradient in g − r, means that a significant number of redder BS stars are likely to scatter into our selection box. Similarly, QSOs have a steep dependence on g − r (see Figure 5), which causes many to scatter into our A-type selection box at fainter magnitudes.

In the following section, our likelihood method is applied only to objects within our selection box. However, by constructing convolved PDF models we are able to take into account the scattering of objects both into and out of this region in ugr space. Our unconvolved models are also defined outside of the ugr selection box (see Equation (7)). Thus, when these models are convolved with photometric uncertainties, the influence of populations which intrinsically lie outside of our ugr selection box (i.e., redder BS stars, and QSOs), are considered in our models. In addition, the convolved models compensate for the affect of fainter stars "leaking" outside of the ugr bounds.

Here, we describe our likelihood analysis used to model the density profile of distant halo stars. We apply our method to SDSS DR9 stars selected in the following magnitude and color range:

$$\begin{eqnarray} 18.5 <& g < 20.5. \nonumber \\ 0.7 <& u-g < 1.6 \nonumber \\ -0.25 <& g-r < -0.1 \end{eqnarray} \tag{ 7 }$$

The number of stars of a particular population (i.e., BHB or BS) in a given increment of magnitude and area on the sky is described by:

$$\begin{eqnarray} \Delta N(m_g\!-\!M_g, \ell, b) &=& \rho (m_g\!-\!M_g,\ell, b) D^3(m_g\!-\!M_g) \nonumber\\ &&\times\frac{1}{5}\mathrm{ln}10 \, \Delta m_g \, \mathrm{cos}b \, \Delta \ell \, \Delta b. \end{eqnarray} \tag{ 8 }$$

Here, we have used Galactic (l, b) coordinates and the heliocentric distance increment ΔD has been converted into the apparent magnitude increment via the relation ΔD = (1/5)ln10 DΔm_g.

We combine Equations (3) and (8) to give the number of A-type stars in a cell of color, magnitude, and longitude and latitude space

$$\begin{eqnarray} \Delta N_{\rm A} &=& p(ugr|{\rm BHB})\Delta N_{\rm BHB}+p({ugr}|{\rm BS})\Delta N_{\rm BS} \nonumber\\ &=& N_{\rm tot}f_{\rm A} \nu _* (ugr, m_g, l,b) \Delta \underline{\mathbf {x,}} \end{eqnarray} \tag{ 9 }$$

where $\Delta \underline{\mathbf {x}}= \mathrm{cos}b \Delta (u-g) \Delta (g-r) \Delta m_g \Delta \ell \Delta b$ is the volume element, and the stellar probability density is

$$\begin{eqnarray} && \nu _* (ugr, m_g, l,b) = \frac{1}{5} \mathrm{ln}10 \, \nonumber\\ && \quad \times \left[f_{\rm BHB} \, p(ugr|{\rm BHB}) \, \rho _{\rm BHB}(m_g\!-\!M_g, \ell, b) D^3_{\rm BHB}(m_g\!-\!M_g)\right. \nonumber\\ && \left.\quad + f_{\rm BS} \, p(ugr|{\rm BS}) \, \rho _{\rm BS}(m_g\!-\!M_g, \ell, b) D^3_{\rm BS}(m_g\!-\!M_g) \right]. \end{eqnarray} \tag{ 10 }$$

Here, the absolute magnitudes of the BHB and BS populations depend on g − r color (M_g = M_g(g − r); see Equation (4)) and f_BHB and f_BS are constants used to ensure that the total number of BHB and BS stars equals the total number of A-type stars (N_A = N_totf_A). We take into account the uncertainty in the BS absolute magnitudes by convolving the number density with a Gaussian magnitude distribution. This distribution is centered on the estimated absolute magnitude ($M^{\rm BS}_g=M^{\rm BS}_g(g-r)$) and has a standard deviation of $\sigma _{M_g}=0.5$.

In our analysis we assume the objects are A-type BHB and BS stars or QSOs. N_tot = N_A + N_Q where N_A = f_AN_tot = N_BHB + N_BS and N_Q = f_QN_tot.

The number of QSOs in a cell of color, magnitude, and longitude and latitude space is

$$\begin{eqnarray} \Delta N_{\rm Q} &=& N_{\rm tot} f_{\rm Q}\nu _{\rm Q} (ugr, m_g, \ell, b) \Delta \underline{\mathbf {x}}, \end{eqnarray} \tag{ 11 }$$

where the QSO probability density, ν_Q is given by Equation (6).

The combined PDF for A-type stars and QSOs is then:

$$\begin{eqnarray} \Delta N_{\rm tot}&=&\Delta N_{\rm A}+\Delta N_{\rm Q}. \nonumber\\ &=&N_{\rm tot}\left((1-f_{\rm Q})\nu _*+f_{\rm Q}\nu _Q\right) \Delta \underline{\mathbf {x}}. \end{eqnarray} \tag{ 12 }$$

Here, we have defined the unconvolved number densities for A-type stars and QSOs. These PDFs are convolved with the ugr, m_g error distributions to give the convolved probability density distribution:

$$\begin{eqnarray} \tilde{\nu }&=& \nu \ast G = \int \int \int \nu (ugr^*, m^*_{g}) \nonumber \\ && \times G(ugr-ugr^*, m_g-m^*_g) {d}(u-g) \, {d}(g-r) \, {d}m_g, \quad\quad \end{eqnarray} \tag{ 13 }$$

where, G(ugr, m_g) is a three-dimensional normal distribution in u − g, g − r, and m_g. The convolved densities are normalized in ugr, m_g, ℓ, and b space over the color and magnitude ranges specified in Equation (7), and over the area of the SDSS DR9 footprint.

The log-likelihood function can then be constructed from the convolved probability density distribution,

$$\begin{eqnarray} \mathrm{log}\mathcal {L} &=& \sum _{i=1}^{N_{\rm tot}} \mathrm{log} \, [ \lbrace (1-f_{\rm Q} ) \tilde{\nu }_* (ugr_i, m_{g,i}, \ell _i,b_i ) \nonumber \\ && +f_{\rm Q}\tilde{\nu }_{\rm Q} (ugr_i, m_{g,i}, \ell _i,b_i ) \rbrace \mathrm{cos}b_i ]. \end{eqnarray} \tag{ 14 }$$

The overall fraction of QSOs, f_Q = 1 − f_A, and relative fraction of BHB stars, f_BHB, are computed iteratively for each set of model parameters from the posterior PDFs:

$$\begin{eqnarray} &&p(\mathrm{QSO}|ugr,m_g)=\frac{N_{\rm Q} \tilde{\nu }_Q}{N_{\rm Q} \tilde{\nu }_Q +N_{\rm A} \tilde{\nu }_*} \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&p(\mathrm{BHB}|ugr,m_g)=\frac{N_{\rm BHB} \tilde{\nu }_{\rm BHB}}{N_{\rm BHB} \tilde{\nu }_{\rm BHB} +N_{\rm BS} \tilde{\nu }_{\rm BS}}. \end{eqnarray} \tag{ 16 }$$

These fractions give the relative contributions of BHBs, BSs, and QSOs in our color–color, magnitude selection box, and ensure the contributions sum to give the total number of stars used in the modeling.

In the following section we outline our model assumptions for the stellar halo density profile. The inner stellar halo density profile (r ≲ 40 kpc) is chosen based on constraints in the literature. We construct a marginal likelihood function by integrating over the adopted range of inner density profile parameters (α₁, α₂, and r_c, see Equations (17) and (18)). The maximum likelihood parameters for the outer stellar halo profile (r_b and α_out, see Equation (17)) are found using a brute-force grid search.

4.1. Model Assumptions

The likelihood method described above is general, and can be applied to any number of model density profiles. Here, we outline the model assumptions applied in our analysis.

The aim of this study is to quantify the outer stellar halo density fall-off. We assume BHB and BS stars follow the same density distribution, which we parameterize as a spherical triple power-law:

$$\begin{equation} \rho (r)_{\rm BHB, BS} \propto \left\lbrace \begin{array}{@{}l@{\quad }l@{}}r^{-\alpha _1} & r \le r_c \\ r^{-\alpha _2} & r_c < r \le r_b \\ r^{-\alpha _{\rm out}} & r > r_b. \\ \end{array}\right. \end{equation} \tag{ 17 }$$

A toy model of this power-law profile is shown in Figure 7 for illustration. In this work, we only consider spherical radial profiles. It is well-known that the stellar halo is flattened in the inner regions (with minor-to-major axis ratio q ∼ 0.6–0.8; e.g., DBE11; Sesar et al. 2011), but it is unlikely that such a flattened profile can exist to large radii. Here, we concentrate on the radial density fall-off, and defer a study of the variation of shape of the stellar halo with radius to future work. In Appendix B we create mock data sets which have flattened stellar halos in the inner regions, and discuss the implications of assuming sphericity at all radii in our modeling procedure.

Zoom In Zoom Out Reset image size

Previous work has shown that within r ≈ 40–50 kpc, the MW stellar density follows a broken power-law (e.g., Bell et al. 2008; Sesar et al. 2011; DBE11). Based on this past work, we assume the following constraints on r_c, α₁, and α₂:

$$\begin{eqnarray} r_c &\in & [20, 30]\, \mathrm{kpc} \nonumber \\ \alpha _1 &=& 2.5 \nonumber \\ \alpha _2 &\in & [3.5, 5.0]. \end{eqnarray} \tag{ 18 }$$

In our analysis, we marginalize over the inner profile parameters. This assumes flat priors over the parameter space given above. Note that the inner-most power-law slope is kept fixed as this has little affect on the outer-most power-law (α_out). The free parameters in our analysis are thus, r_b and α_out, and we consider values in the range: r_b ∈ [30, 70] kpc and α_out ∈ [2.0, 10.0].

4.2. Tests with Mock Data

To demonstrate the ability of our modeling technique, we apply our likelihood method to "mock" data. For our mock data, we assume a (fixed) inner profile with r_c = 25 kpc, α₁ = 2.5, and α₂ = 4.0, and consider three different outer halo models: a "shallow" model with r_b = 40 kpc, α_out = 3.0, a "constant" fall-off model with α_out = α₂ = 4.0, and a "steep" model with α_out = 6, r_b = 50 kpc. In all mock data sets, we adopt overall population fractions of f_BHB = 0.18 and f_Q = 0.23.

The following steps are applied to generate the mock data.

• 1.  
  A-type stars (BHB and BS) and QSOs are drawn from the (unconvolved) model PDFs defined in Equations (6) and (10) using an acception-rejection algorithm. Objects are generated in ugr, m_g, ℓ, b space from uniform color, magnitude distributions, and l, b are drawn randomly from the surface of a sphere.
• 2.  
  For the A-type stars, g-band magnitudes are converted to heliocentric distances using absolute magnitude–color relations appropriate for each stellar population. The BS absolute magnitudes are scattered about their mean relations assuming $\sigma _{M_g} \sim 0.5$.
• 3.  
  Only high latitude objects (with |b| > 30°), inside of the SDSS DR9 footprint are considered.
• 4.  
  The three-dimensional error distribution in ugr, m_g space appropriate for our SDSS DR9 sample (see Figure 2) is applied to the mock data. After photometric scattering, only objects lying within the bounds defined in Equation (7) are considered.
• 5.  
  Our mock data sets are generated with the same number of stars as our SDSS DR9 sample with known substructures removed (N = 5213, see Section 4.3).

The magnitude distributions of our three models ("shallow": dashed green, "constant": solid blue, "steep": dot-dashed red) are shown in Figure 8. The overall magnitude distributions for the three models are similar, but there are clear differences between the BHB star distributions; this is not surprising given that the distance range of the BS stars generally lie within the (fixed) inner density profile (e.g., approximately 77% of BS stars are within r = 30 kpc).

Zoom In Zoom Out Reset image size

The results of applying our likelihood method to the mock data are shown in Figure 9. In the left-hand panel we show the likelihood contours in α_out, r_b space. The filled and unfilled contours indicate 1σ and 2σ confidence regions respectively. In all cases, our method is able to reproduce, within the uncertainties, the true density profiles. In the right-hand panel we show the maximum likelihood α_out values for different fixed values of break radius r_b. The lines indicate the median values, and the shaded regions encompass the 1σ confidence regions. The "steep" models show a characteristic steepening as the adopted break radius is increased. Such models are clearly distinguished from shallower profiles.

Zoom In Zoom Out Reset image size

In the above exercise we adopt the same inner profile as the input mock data (α₁ = 2.5, α₂ = 4.0, r_b = 25 kpc). However, when we apply our method to the SDSS DR9 data we marginalize over a wide range of inner profile parameters (see Equation (18)). In Appendix B we show that this flexibility can compensate for biases induced by assuming an inaccurate inner stellar halo profile.

4.3. Treatment of Known Substructures

In our analysis, we consider the effect of known structures on our density profile estimates. There are two large structures in the regime of our sample which could affect our results: the Virgo overdensity (Jurić et al. 2008) and the Sagittarius (SGR) stream (Ibata et al. 1995).

The Virgo overdensity is located at high latitudes and mainly affects the brightest BS stars in our sample (in the distance range 10 ≲ D/kpc ≲ 20). We isolate stars belonging to Virgo by applying the following cut in Galactic coordinates (see Bell et al. 2008; DBE11):

$$\begin{equation} 0 < X < 30, \quad \quad X=63.63961\sqrt{2(1 - \sin b)}. \end{equation} \tag{ 19 }$$

Stars belonging to the SGR stream are present over the full magnitude range of our sample (see Figure 13, Belokurov et al. 2014). Thus, we locate possible SGR stars according to their position on the sky (see Deason et al. 2012, Figure 7).

Finally, a small fraction of distant BHB stars in our sample coincide with two known dwarf galaxies; Sextans (D ∼ 90 kpc) and Ursa Minor (D ∼ 60 kpc). Stars in the regions of these dwarfs (N ∼ 200) are excluded.

Our selection of A-type stars gives N_tot = 10, 787 objects including SGR and Virgo, and N_tot = 5213 when objects in the regions of these known overdensities are excluded. Below, we apply our likelihood analysis to our SDSS DR9 sample both with and without these large substructures.

In this section, we apply our likelihood technique to our sample of A-type stars selected from SDSS DR9. In Figure 10 we show the likelihood results. In the left-hand panel, we show the 1σ and 2σ confidence levels for the outer slope (α_out) and break radius (r_b). The dashed red and solid black lines show the results both with and without SGR and Virgo. In the right-hand panel we show the maximum likelihood outer slope for different fixed values of break radius.

Zoom In Zoom Out Reset image size

There is a strikingly steep fall-off (α_out ∼ 6–10) in the stellar halo density beyond ∼50–60 kpc. This holds true even when the SGR and Virgo overdensities are included in the analysis. The implications of this result for the accretion history and mass profile of the Milky Way are discussed in Section 6. We note that the location of a "break" in the stellar density at r ∼ 50–60 kpc is coincident with the apocenter of the SGR leading arm (see Belokurov et al. 2014 Figure 4). Deason et al. (2013a) showed that "breaks" in the stellar halo density are strongly linked to the apocenters of their accreted constituents; thus, the location of a break close to the apocenter of SGR agrees very well with this hypothesis.

In our analysis, the fraction of QSO contamination (f_Q) and the overall BHB star fraction (f_BHB) are found iteratively for each model PDF. In Figure 11 we show the variation of these fractions with our model parameters. Our maximum likelihood model parameters are summarized in Table 1. Here, we also give the maximum likelihood parameters over the four-dimensional space (r_c, α₂, r_b, α_out). The maximum likelihood inner profile parameters (r_c, α₂) listed here, dominate over the likelihood when we marginalize over a range of parameters (see Equation (18)).

Zoom In Zoom Out Reset image size

Table 1. Maximum Likelihood Results

rb (kpc)	αout	fQ	fBHB	$\Delta \mathrm{log}\mathcal {L}$a
Exc. Virgo and SGR
30	$4.8^{+0.3}_{-0.2}$	0.234	0.193	−4.9
35	$4.8^{+0.2}_{-0.3}$	0.233	0.192	−4.5
40	$4.6^{+0.3}_{-0.2}$	0.231	0.197	−4.7
45	$4.8^{+0.5}_{-0.3}$	0.232	0.194	−5.3
50	$6.0^{+0.6}_{-0.9}$	0.235	0.180	−4.6
55	$7.0^{+0.8}_{-0.8}$	0.235	0.178	−2.4
60	$8.4^{+1.2}_{-0.8}$	0.236	0.178	−0.6
65	>7.1b	0.236	0.181	0.0
70	>7.1b	0.235	0.189	−0.3
(rb, αout)c$=(\, 65^{+5}_{-6} \, \mathrm{kpc}$,  >6.2b)
(rc, α1, α2, rb, αout)MLd =( 25 kpc, 2.5, 4.5, 65 kpc, 10 )
Inc. Virgo and SGR
30	$4.4^{+0.2}_{-0.2}$	0.198	0.249	−22.2
35	$4.4^{+0.2}_{-0.3}$	0.195	0.248	−22.8
40	$4.6^{+0.2}_{-0.2}$	0.194	0.245	−20.6
45	$5.4^{+0.2}_{-0.3}$	0.196	0.230	−8.2
50	$6.0^{+0.5}_{-0.3}$	0.195	0.228	−2.2
55	$7.2^{+0.5}_{-0.5}$	0.196	0.223	0.0
60	$9.0^{+0.8}_{-0.8}$	0.196	0.221	−2.6
65	$7.8^{+1.2}_{-0.8}$	0.194	0.242	−10.9
70	$8.4^{+1.0}_{-1.0}$	0.193	0.248	−14.3
( rb, αout)c$=(55^{+4}_{-3} \, \mathrm{kpc}, 7.2^{+1.6}_{-0.7}$ )
(rc, α1, α2, rb, αout)MLd =( 30 kpc, 2.5, 3.5, 55 kpc, 7.2 )

Notes. ^aDifference in log-likelihood from maximum likelihood value. ^b2σ lower limits. ^cJoint maximum likelihood result for r_b and α_out after marginalizing over the inner stellar halo parameters r_c and α₂. ^dThe maximum likelihood parameters from the four-dimensional space of r_c, α₂, r_b, α_out. Note, α₁ is kept fixed (see Equation (18)).

Download table as:  ASCIITypeset image

In Figure 12 we show the g-band magnitude distribution of the data and best-fit model. The blue/red shaded regions show the best-fit models when SGR and Virgo are excluded/included respectively. In the left panel we show all the stars in our ugr selection box. In the right panel we only show stars with u − g > 0.9 to reduce the influence of QSOs on the g-band magnitude distribution. There is good agreement between the data and models, especially when known overdensities are excluded. Finally, we show in Figure 13 the residuals of our models and data on the sky (in equatorial coordinates). We split into two magnitude bins; "bright" (18.5 < g < 19.5) and "faint" (19.5 < g < 20.5). The increasing dominance of SGR is evident in these two panels. However, we note that away from the SGR stream, the residuals are close to zero.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.1. Milky Way Accretion History

Our finding of a strikingly sharp drop in stellar density beyond r ∼ 50–60 kpc may have important implications for the accretion history of the stellar halo. In particular, our results suggest that, other than the relatively recent accretion of the SGR dwarf, the "cannibalistic" past of the Milky Way likely subsided several gigayears ago.

To illustrate the dependence of the outer stellar halo slope on its past accretion history, we compare with the Bullock & Johnston (2005) stellar halo models. This suite of 11 simulated stellar halos are built up entirely from the disruption of dwarf galaxies. The accretion history of each M_vir ∼ 1.4 × 10¹² M_☉ mass halo is generated at random using semi-analytic merger trees appropriate for a ΛCDM cosmology. In the left panel of Figure 14 we show the density profiles of the 11 halos between 50 < r/kpc < 100. It is clear that there is no "universal" outer halo fall-off, and there is a wide variation in the density profiles. In the right panel of this figure we show the outer stellar density slope against the average time at which the stars in the radial regime 50 < r/kpc < 100 became unbound from their parent dwarf (T_ub). The filled gray region indicates the approximate slope for the Milky Way, α_out[r_b = 50 kpc] = 6.

Zoom In Zoom Out Reset image size

Despite some scatter, it is clear that halos with shallower slopes tend to have more recent accretion activity. It is worth noting that we have made no attempt to "excise" substructure from these simulated halos, and this will likely lead to steeper profiles in some cases. However, the general trend indicates that beyond r ∼ 50 kpc in the Milky Way halo, the "field" halo stars were likely stripped from dwarfs that were accreted a long time ago (>6 Gyr).

Finally, we note that the Eris simulation (Guedes et al. 2011), one of the highest resolution hydrodynamical simulations of the formation of a M = 8 × 10¹¹ M_☉ late-type spiral, also has a very steep fall-off in stellar density beyond r ∼ 60–70 kpc (see Rashkov et al. 2013 Figure 2). This simulation, which has been successful in matching several Milky Way properties, has an early accretion history and high concentration (c_vir ∼ 24). This adds further weight to our deductions from the Bullock & Johnston (2005) simulations, that the Milky Way halo has undergone a relatively quiescent accretion history over the past several gigayears.

6.2. Milky Way Mass

6.3. Comparison with M31

We model the density distribution of distant BHB and BS halo stars using SDSS DR9 photometry, with the aim of measuring the outer slope of the Milky Way stellar halo density profile beyond r ∼ 50 kpc. We construct number density PDFs in ugr, m_g space, and include contributions from QSO contaminants. Our PDF is convolved with the ugr, m_g error distribution to take into account the significant photometric uncertainties at faint magnitudes. We fix the QSO number density using the QSO model developed by Bovy et al. (2011a), and allow the stellar halo profile within r ∼ 40–50 kpc to lie within an observationally motivated parameter space. The outer halo model parameters are identified by modeling the stellar distribution in u − g, g − r, m_g, ℓ, b space. We test our method on simulated catalogs of BHBs, BSs, and QSOs, and demonstrate that the properties of the distant halo can be recovered with sufficient accuracy.

We apply our likelihood analysis to high latitude (|b| > 30 deg) SDSS DR9 stars in the color and magnitude range; 0.7 < u − g < 1.6, −0.25 < g − r < −0.1, and 18.5 < g < 20.5. With this selection, BHB and BS stars span a heliocentric distance range: 10 ≲ D_BS/kpc ≲ 75, 40 ≲ D_BHB/kpc ≲ 100. We identify stars coincident on the sky with the known substructures Virgo and SGR, and apply our analysis both including and excluding these stars. Our analysis assumes (1) stellar halo sphericity at large radii, (2) an inner stellar halo (r ≲ 40 kpc) density parameterization consistent with current constraints in the literature, and (3) BHB and BS intrinsic color distributions that remain the same throughout the halo.

The relative contributions of A-type stars (BHB and BS) and QSOs are computed iteratively from the convolved PDFs for each set of model parameters. In our selection box 0.7 < u − g < 1.6, −0.25 < g − r < −0.1, and magnitude range 18.5 < g < 20.5, we find a QSO contamination fraction of f_Q ∼ 0.2 and a BHB fraction of f_BHB ∼ 0.2; these fractions have a weak dependence on our model parameters.

After excluding known substructures, we find that very steep outer halo profiles are preferred, with α_out ∼ 6 beyond r = 50 kpc. Even when SGR and Virgo stars are included in the analysis, we find very steep outer profiles. There is evidence for a break in the stellar density at r_b ∼ 50–60 kpc, which is coincident with the apocenter of the SGR leading arm.

We compare our results to the predictions of simulated stellar halos. The Bullock & Johnston (2005) suite of halos, built up purely from the accretion of dwarf galaxies, have outer profile slopes which depend on the accretion history of the halo; steeper outer slopes suggest earlier accretion epochs than shallow slopes. Thus, our finding of a very steep outer halo profile argues that, apart from the relatively recent accretion of SGR, the majority of the Milky Way stellar halo was built up from relatively early accretion events (T > 6 Gyr ago). This is in contrast to the M31 stellar halo which has a much shallower density slope out to r ∼ 200 kpc (α ∼ 3–3.5; Gilbert et al. 2012; Ibata et al. 2014), and thus presumably has had a more active late-time accretion history.

The density profile of the Milky Way stellar halo is an important ingredient for dynamical mass estimates of the Galaxy. Until now, the unknown stellar density slope beyond ∼50 kpc has proved to be a troublesome bottleneck in constraining the total mass out to large distances. Our finding of a very steep outer halo slope may have important implications for studies utilizing the kinematics of halo stars to estimate the total mass of the Milky Way. The measurement we report here, in combination with constraints on the halo star velocity anisotropy from upcoming surveys (such as Gaia) will, undoubtedly, significantly reduce the uncertainty surrounding dynamical mass measurements of the Milky Way.

We thank Jay Farihi for his help with the WD models. A.J.D. is currently supported by NASA through Hubble Fellowship grant HST-HF-51302.01, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. V.B. thanks the Royal Society and the European Research Council for financial support. We thank the Aspen Center for Physics and the NSF grant No. 1066293 for hospitality during the conception of this paper. We also thank an anonymous referee for useful comments.

We select WDs by applying the XDQSO algorithm (see Section 3.1.1) to our high latitude (|b| > 30°) SDSS DR9 photometry in the magnitude range 18.5 < g < 20.5, and only consider blue stars (u − g < 0.7) with a low QSO probability (p_qso < 0.1).

The left three panels of Figure 15 show the distribution of WD stars in u − g color space for three g − r color bins. Two populations are evident: these are the H-dominated (redder u − g) and He-dominated (bluer u − g) populations, which we refer to as DA and DB type respectively. The solid lines show a double Gaussian fit to these populations. The vertical dashed (dot-dashed) lines indicate the loci of model WDs with surface gravity log(g_s) = 8.0(7.5). These model WDs derive from the Montreal WD atmosphere group⁸ (Holberg & Bergeron 2006; Kowalski & Saumon 2006; Tremblay et al. 2011; Bergeron et al. 2011), who provide synthetic model sequences for WDs with SDSS ugriz photometry.

Zoom In Zoom Out Reset image size

In the two right panels of Figure 15 we show the ridgelines of the WD populations in ugr space based on the double Gaussian fits. The dashed lines indicate the ridgelines predicted by models with log(g_s) = 7.5–8.0. The agreement is good, especially for the DA-type dwarfs which dominate the WD population. We define the WD ridgelines with polynomials (see Equation (1)):

$$\begin{eqnarray} (u-g)^{0}_{\rm DA}&=& 0.400+0.046(g-r)+1.186(g-r)^2 \nonumber\\ && -2.227(g-r)^3-42.785(g-r)^4 \nonumber \\ (u-g)^0_{\rm DB}&=& 0.245+1.098(g-r)-0.412(g-r)^2 \quad \end{eqnarray} \tag{ A1 }$$

for −0.25 < g − r < −0.1.

We calculate the intrinsic spread of the two populations about their ridgelines to be σ_{DA, 0}(u − g) = 0.060 and σ_{DB, 0}(u − g) = 0.075. In a similar fashion to the A-type stars, we assume Gaussian distributions about the ridgelines:

$$\begin{eqnarray} && p(u \,{-}\, g|{\rm DA, {\rm {g}} \,{-}\, {\rm {r}}})\propto \mathrm{exp}\left(\!{-}\frac{\left[(u \,{-}\, g) \,{-}\, (u \,{-}\, g)_{\rm DA}^0\right]^2}{2\sigma _{\rm DA}^2} \!\right), \nonumber \\ && p(u \,{-}\, g|{\rm DB, {\rm {g}} \,{-}\, {\rm {r}}})\propto \mathrm{exp}\left(\!{-}\frac{\left[(u \,{-}\, g)\,{-}\,(u \,{-}\, g)_{\rm DB}^0\right]^2}{2\sigma _{\rm DB}^2} \!\right). \quad\quad \end{eqnarray} \tag{ A2 }$$

We assume a constant intrinsic g − r distribution, so the color-based probabilities of class membership are then: p(ugr|DA)∝p(u − g|DA, g − r), p(ugr|DB)∝p(u − g|DB, g − r). We fix the fraction of DA-type WDs (assuming just DA and DB types) to be f_DA = 0.7.

In the right-most panel of Figure 15 we show the absolute magnitude–color relation for the model WDs (with log(g_s) = 7.5–8.0). The DA and DB-types have similar absolute magnitudes. We adopt the weighted mean of these relations (f_DA = 0.7, f_DB = 0.3) as the average WD absolute magnitude relation, and assume a 0.5 mag spread to account for uncertainties in log(g_s) (e.g., log(g_s) = 7.5(8.0) models have brighter(fainter) absolute magnitudes by ∼0.5 dex):

$$\begin{equation} M_{g(\rm WD)}= 12.249 + 5.101(g-r), \end{equation} \tag{ A3 }$$

where, $\sigma _{M_g(\rm WD)} \sim 0.5$.

Finally, we fix the WD density distribution assuming a disk distribution of stars. We use the disk density profile found by Jurić et al. (2008), which assumes an exponential profile and has contributions from thin and thick disk populations:

$$\begin{eqnarray} \rho _{\rm thin} &=& \mathrm{exp}\left(R_0/L1\right) \mathrm{exp}\left(-R_{\rm WD}/L1-|z_{\rm WD}+z_0|/H1\right) \nonumber \\ \rho _{\rm thick} &=& \mathrm{exp}\left(R_0/L2\right) \mathrm{exp}\left(-R_{\rm WD}/L2-|z_{\rm WD}+z_0|/H2\right) \nonumber \\ \rho _{\rm WD} &=& \rho _{\rm thin}+0.12\rho _{\rm thick}, \end{eqnarray} \tag{ A4 }$$

where, H1 = 0.3, L1 = 2.6, H2 = 0.9, L2 = 3.6, z₀ = 0.025 kpc, R₀ = 8.5 kpc.

Figure 4 in the main text shows that our WD models predict a very small fraction of WDs in our A-type star selection box, so we do not consider their contribution in our modeling procedure.

Our analysis assumes spherical stellar halo density profiles. Here, we consider the implications of this assumption for non-spherical profiles. Mock data is generated, as described in Section 4.2, but our model stellar halo profiles are given a minor-to-major axis ratio, q, which varies with radius. We adopt the following parameterization for q:

$$\begin{equation} q(r)=q_0 \sqrt{\frac{r^2+r_s^2}{\left(q_0 r\right)^2 +r_s^2}} \end{equation} \tag{ B1 }$$

and set q₀ = 0.6 and r_s = 15 kpc. The minor-to-major axis parameter thus varies smoothly from q ∼ 0.6 at small radii to q ∼ 1.0 at large radii (see inset in bottom right panel of Figure 16). Three mock data sets are generated with the same parameters as described in Section 4.2. We apply our likelihood analysis to this simulated data assuming sphericity at all radii.

Zoom In Zoom Out Reset image size

The results of this exercise are summarized in Figure 16 (cf. Figure 9). The "constant" (blue), "steep" (red), and "shallow" (green) toy models are shown in the top, middle, and bottom rows, respectively. Solid filled regions and solid lines indicate the 1σ and 2σ confidence regions when a spherical model (q = 1) is used in the likelihood analysis with inner stellar halo parameters: α₁ = 2.5, α₂ = 4.0, r_c = 25 kpc. In all cases, the resulting outer stellar density parameters (r_b and α_out) are biased toward shallower profiles: typically α_out is 0.5 dex too shallow. Our flattened model mainly affects the BS stars at small radii; when forced to fit to a spherical model, the BS distribution appears shallower than the input inner density profile (α₂(input) = 4.0). The best-fit model, with only α_out and r_b as free parameters, compensates for this bias by making α_out slightly shallower.

The line-filled regions and dashed lines indicate the 1σ and 2σ confidence regions when a spherical model (q = 1) is used in the likelihood analysis with inner stellar halo parameters: α₁ = 2.5, α₂ = 3.5, r_c = 25 kpc. Thus, this model adopts a shallower inner profile (α₂ = 3.5) than the case described above. In this case, we are able to reproduce the correct outer stellar halo parameters with reasonable accuracy. By using a shallower inner profile in the modeling, we have "compensated" for the affect of flattening. So, we are able to reproduce the correct stellar density profile at large radii, even though we have neglected the affects of flattening.

When we apply our analysis to the SDSS DR9 data, we marginalize our likelihood distribution over a wide range of inner profile parameters (see Equation (18)). This ensures that our inner profile is flexible enough to compensate for affects such as a flattened inner profile. The dotted black lines in the left panels of Figure 16 indicate the 1σ contours after marginalizing over the two inner density profile models. It is clear that the marginalized likelihood is dominated by the (higher likelihood) model which is able to reproduce the correct outer stellar density slope. This gives us confidence that the resulting outer stellar density parameters are robust to variations in the inner stellar halo profile.

Note that we also use mock data to test how variations in α₁ may affect our results (this is kept fixed at α₁ = 2.5 in our analysis). We find that the outer profile parameters are less sensitive to variations in α₁ than α₂ or r_c, and we find little difference if α₁ is changed by ∼ ± 0.5 dex.