ON THE SHOULDERS OF GIANTS: PROPERTIES OF THE STELLAR HALO AND THE MILKY WAY MASS DISTRIBUTION

We mine the ninth SDSS data release DR9 (Ahn et al. 2012) to construct our giant star catalog. The first set of selection criteria we impose to prepare the catalog are

where m_g and m_r are the extinction-corrected magnitudes. The metallicity [Fe/H] and the stellar parameters we use, i.e., surface-gravity log g and effective temperature T_eff, are the ones labeled in the SDSS DR9 as "ADOP," meaning the average of various estimators that the SSPP ³ pipeline uses, and also are the ones recommended by Yanny et al. (2000). The quadratic color-metallicity cut in Equation (3), devised by Xue et al. (2014) and inferred from An et al. (2008), is basically meant to eliminate the contamination from red horizontal branch and red clump stars. Note that we take a slightly conservative cut on log g to minimize the contamination of dwarfs. Additionally, we apply a second set of selection criteria given by

as a quality control cut. This is to ensure the accuracy and the reliability of the stellar parameters and radial velocity values that SSPP provides. Combination of the above given sets of selection criteria yields 5330 candidate giants. For further diagnostics about selection criteria and candidate giants, see Appendix A.

As shown in the next section, we estimate the most probable distance of the stars, which gives the height, z, of the stars from the Galactic plane. We then impose the condition |z| > 4 kpc to select the halo stars. This leaves us with 5140 stars.

Correct distance measurements of the giant stars are critical for studying the kinematics of the halo and also for modeling the mass. For a correct treatment of the observational errors, we set up the distance estimation in a Bayesian framework; the calculations are given in full in Appendix B. The procedure we follow is the same as given by Xue et al. (2014). The essence of the exercise is that for each star with some set of observables, say, S = {m, c, [Fe/H]}, we obtain a corresponding absolute magnitude by matching it to color-metallicity fiducials of red giant sequences of clusters. In the end, with the inferred absolute magnitude and a given apparent magnitude, we use the standard photometric parallax relation to compute the distance for a star. Instead of a single valued number, the setup allows us to compute a full probability distribution or posterior distribution p(μ|S) of the distance modulus μ.

Studies of the morphology of the Milky Way halo suggest that a good fit to the stellar halo density distribution is a double power law with a shallow slope inside a break radius (r_b) and a sharp fall-off outside. For example, analyzing the main-sequence turnoff stars, Bell et al. (2008) conclude that a power law, ρ∝r^−α, with index of α = 2 for the inner region and 4 in the outer region with a break at r ≲ 20 kpc is a reasonable representation. Similar conclusions were also made by Watkins et al. (2009) in their studies of RR Lyrae stars out to 100 kpc and by Deason et al. (2011) for BHB stars out to 40 kpc. Apart from SDSS, the study by Sesar et al. (2011) of the main-sequence turnoff stars obtained from the Canada–France–Hawaii Telescope Legacy Survey suggests a slightly shallower fall of 3.8 beyond the break radius r = 28 kpc. More recently, Sesar et al. (2013) using variable stars suggest a broken power law but with a much smaller break radius of ∼16 kpc. In agreement with Watkins et al. (2009), Deason et al. (2011), here we adopt a double power-law fit to the halo density with an inner slope of 2.4 and an outer slope of 4.5 with break at radius r_b. As we will discuss later, the drop seen in the σ_r profile at r ≈ 20 kpc seems likely to be a consequence of the break in the density. Hence, to further investigate this issue, we keep the break radius (r_b) as a free parameter of our model.

The outermost (r ≳ 100 kpc) region of the halo is diffuse and highly structured (Sesar et al. 2007; Sharma et al. 2011b). Only a handful of stars have been observed at such large distances, and so the density profile is largely unknown. With a catalog of a comparatively larger number of outermost halo stars, we investigate what the decline of σ_r(r) after r ≳ 90 kpc tells us about the density profile. For this, we investigate a truncated model of the outer halo with a truncation radius r_t. A similar approximation has also been made by Dehnen et al. (2006), but with a forced sharp truncation at r = 160 kpc. However, we soften the truncation by using an exponential functional form that has a tunable parameter Δ that determines the strength of the fall (see also Sharma et al. 2012). We then let the MCMC likelihood fit determine both the position of r_t and the scale length of the fall Δ. The first row in Table 1 shows the form of the density profile that we adopt. Note that we assume the logarithmic slope at r = r_t to be continuous. This gives = (r_t/Δ) − 4.5, and the slope at r ⩾ r_t as − r/Δ.

The velocity anisotropy (β) is the most uncertain quantity that enters the Jeans equation, and only recently has β(r) been measured directly. The full phase-space studies of the subdwarfs in Smith et al. (2009a) and main-sequence stars in Brown et al. (2010) suggest that the halo within d < 10 kpc is radial with β ≈ 0.7. Recently, K12 studied the line-of-sight velocity of BHB stars and found that at a similar distance range (r < 12 kpc) the halo is radial β ∼ 0.5. Moreover, K12 also found that the β(r) profile of the stellar halo has features that include a tangential dip (β = −1.2) at r ∼ 17 kpc. In the range 19 ≲ r/kpc ≲ 25, K12 find the β to be consistent with zero given the uncertainty. Further confirmation of this comes from Deason et al. (2013), who used the HST proper-motion information of the main-sequence stars to find that at similar distance r ∼ 25 kpc the halo is isotropic $\beta =0.0^{+0.4}_{-0.2}$. For r > 25 kpc, so far there has been no direct measurement of β(r). Therefore, in this regime we have no choice but to assume some model for β. A trend that β(r) rises from 0 and attains a positive value at large distances has been reported (see Figure 13 of K12) for the ΛCDM stellar halos of Bullock & Johnston (2005) and also for the halo stars in the cosmological simulations of Sales et al. (2007). We assume a similar trend here, i.e., at distance r > r₂, β(r) is constant and equal to the β_s value. However, in between the maximum observed radius r = 25 kpc and the distance r₂, β(r) is assumed to be linear. This is more physical than introducing an abrupt transition like a step function. In other words, r₂ determines the slope of the β(r) profile joining β_s and β|_{r = 25 kpc}. Here we assume r₂ to be 50 kpc. This transition is also consistent with the results from the simulation. However, in both of the above-mentioned simulations the transitions cease at a much closer radius than we assume here. The formula shown in the second row of Table 1 summarizes our model for the anisotropy profile.

The next quantity required for the σ_r modeling is the potential, Φ(r). We essentially construct a model of the Galactic potential with three components: an oblate spheroidal bulge, an axisymmetric disk described by Miyamoto & Nagai (1975), and a spherical dark halo described by a ΛCDM motivated Navarro–Frenk–White (NFW; Navarro et al. 1996) profile. The formulae are given in Table 1.

The Galactic dark matter halo assumed to follow an NFW profile is characterized by two parameters, the virial mass M_vir and the concentration c (third row of Table 1). Here the overdensity of the dark matter compared to the average matter density, δ_th, is considered to be 340 (Bryan & Norman 1998). The values for Hubble constant H₀ = 70.4 kms⁻¹ Mpc⁻¹ and the matter density of the universe Ω_m = 0.3 are taken from Komatsu et al. (2011); ⁴ in other words, we assume that the mean density is Ω_mδ_th times the cosmological critical density.

The Galactic disk is thought to have two major components, i.e., a thick and a thin disk (e.g., Gilmore & Reid 1983). Generally, the disk is modeled as exponential in a sense that the surface density falls exponentially as a function of distance from the center of the Galaxy R and height from the Galactic midplane z. Since here we use the spherical Jeans equation, we consider an analytic and easy-to-use three-dimensional model of the disk, which is provided by a flattened disk of Miyamoto & Nagai (1975) type. Its functional form is given in the second-to-last row of Table 1. There, a and b are its scale lengths and M_disk is the mass.

The Galactic bulge is assumed to be spheroidal. This provides a reasonable axisymmetric approximation for the bulge seen from COBE/DIRBE near-infrared data. It is also similar to the axisymmetric approximation of the Bissantz & Gerhard (2002) model considered in McMillan (2011). The last row of Table 1 presents the mass density of the assumed model for the bulge and is taken from Equations (2.207a) and (2.207b) in Binney & Tremaine (2008). To compute the force generated by such a spheroidal system, we use Equations (2.129a) and (2.129b) of Binney & Tremaine (2008). The values for some of the bulge parameters are kept fixed, such as oblateness parameter q_b and power-law index α_b. These were adopted from Section 10.2.1 of Binney & Merrifield (1998). We found that at a distance greater than four times the scale length s, the spheroidal bulge contribution to the overall potential is similar to that of a point mass.

Both the disk and bulge models are functions of radius and polar angle. But here we work in a framework of the spherical Jeans equation. Hence, we only consider the radial component of the force due to disk and bulge, i.e., along the basis vector $\widehat{e_r}$, which we average over the spherical shells. The average force exerted on a unit mass at given radius r due to all three components of the Galaxy is modeled by

The dominant contributor to the overall potential of the Galaxy in the innermost region r ≲ 5 kpc is the bulge, in 5 kpc < r < 15 kpc is the disk, and at even larger distances is the halo. Therefore, fixing a potential of any component would systematically alter the contribution of the others. Thus, we keep the parameters free and allow the data to resolve the degeneracies.

Before proceeding with the fit, we first study the effect of density and anisotropy profiles in the σ_r(r) model. The model σ_r(r) is obtained by substituting the above-described density, anisotropy, and potential profiles in Equation (9). In Figure 3(a), we demonstrate the role of the assumed density profile for the case of β = 0.5, c = 10, and M_vir = 10¹² M_☉. The solid (dotted) line is the case of a single power law with a slope of 2.4 (4.5), whereas the dashed line is the case of a double power-law density profile with an inner slope 2.4 and an outer slope 4.5, with the break being at r_b = 50 kpc. In the figure the σ_r for the double power-law case attains the values for single power-law cases in the inner and outer parts with a sharp transition ceasing exactly at the break radius r_b. Furthermore, the dashed-dotted line in the figure is the same as in the case of the dashed line but with an exponential fall-off beyond r_t = 120 kpc with parameter Δ = 1.1 kpc. Note that the σ_r for the dashed-dotted line and dashed line cases are the same out to r_b = 50 kpc, but beyond that the dashed-dotted line declines further owing to added exponential fall-off in the density profile. Since the slope in the latter case is a function of r, the transition is smoother than one near the first break r_b. For r > r_t, σ_r decreases gently. This suggests that the break in the σ_r(r) profile is a result of the break, r_b, in the density distribution.

Zoom In Zoom Out Reset image size

In Figure 3(b), we demonstrate the effect of the underlying β profile for the case of c = 10, M_vir = 10¹² M_☉, density power-law index of 2.4 (inner region) and 4.5 (outer region), and r_b = 23 kpc. The β(r) plotted here are taken from the second row of Table 1. The dotted line and dashed-dotted lines show β(r) with β_s = 0.5 and r₂ = 75 kpc and 50 kpc, respectively. These two runs can be compared to see the effect of r₂ on the σ_r(r) model. Clearly, the chosen value of r₂ that determines the distance at which β saturates and attains constant β_s value is just the measure of the slope of the transition. Also, as expected, in the case of r₂ = 75 kpc the slope is smaller than in the case of r₂ = 50 kpc as it attains the β_s at larger r. The effect of r₂ seems minimal in the overall σ_r(r) profile. Hence, we keep it fixed at 50 kpc for the rest of the analysis. Furthermore, the above two runs can be compared with the dashed line to see the effect of adopted β_s, which is 0.5 in the former runs and 0.9 in the latter one. As expected, larger β_s = 0.9 means that the σ_r(r) shifts up and vice versa. The β_s can systematically shift the σ_r(r) profile and bias the mass estimate; we keep it free.

The model parameters we are interested in measuring are the stellar halo density power-law break radius (r_b), truncation radius (r_t), truncation softening (Δ), and maximum anisotropy (β_s); dark matter halo concentration (c) and the virial mass (M_vir); disk scale lengths (a and b) and mass (M_disk); and bulge scale length (s) and mass (M_bulge). We consider flat priors for all the parameters in the following range: r_b ∈ [8, 30] kpc, r_t ∈ [60, 140] kpc, Δ ∈ [0, 20] kpc, β_s ∈ [ − 5, 1], c ∈ [1, 60], M_vir ∈ [0.05, 3] × 10¹² M_☉, M_disk ∈ [0.1, 5] × 10¹¹ M_☉, a ∈ [0.1, 12] kpc, b ∈ [0.01, 0.5] kpc, M_bulge ∈ [0.1, 3] × 10¹⁰ M_☉, and s ∈ [1, 3] kpc.

Next, we use the MCMC and likelihood maximization to compute the model parameters of interest. The likelihood function we use is

where θ = {r_b, r_t, Δ, β_s, c, M_vir, M_disk, a, b, M_bulge, s} is the set of model parameters we explore. The last term is the prior on ω_LSR. The first term is

In this model σ_r(r_k, θ) is given by Equation (9). The probability p(σ_r|r_k) is the posterior distribution of σ_r parameter at the kth node and is obtained from Equation (8) in Section 3, in other words, the distribution of each data point in Figure 1. This setup avoids assuming gaussian errors for the measured values; instead, we already have the full probability distribution of σ_r at each r, and we make use of this.

Our catalog only contains the halo stars, and with the halo stars alone we cannot construct the rotation curve for the inner region of the Galaxy. Fortunately, the shape of the rotation curve or the circular velocity v_circ in the innermost region of the Galaxy r < R_☉ kpc, where the bulge and the disk dominate, can be computed from alternative measures such as using tangent point velocities (e.g., Malhotra 1995; McClure-Griffiths & Dickey 2007) or the gas rotation curve (e.g., Sofue et al. 2009). Here we use the terminal velocity curves from Malhotra (1994, 1995), as done, e.g., in Dehnen & Binney (1998), Widrow et al. (2008), McMillan (2011).

By measuring the terminal velocity along all lines of sight between Galactic longitudes |l| < π/2 for latitude b = 0°, it is possible to derive a measurement of the rotation curve of the inner Galaxy. Assuming that the Galaxy is axisymmetric, v_term as a function of l is given by

where v_LSR = v_circ(R_☉). Now we define our second term in likelihood as

There will be effects of nonaxisymmetry of the Galaxy and noncircular motion of the interstellar medium on the $v_{\rm term} ^{\rm data}$. To take into account these effects, following Dehnen & Binney (1998), we assume $\sigma _{v_{\rm term} } = 7 \ {\rm km\ s^{-1}}$ and avoid the region affected by the bar by only using data with |sin l| > 0.3. The data with assumed uncertainties are shown with black points in Figure 4.

Zoom In Zoom Out Reset image size

An additional prior we impose is on v_LSR. There is a wide variation in claims about v_LSR at R_☉ ranging between 184 km s⁻¹ (Olling & Merrifield 1998) and 272 km s⁻¹ (Méndez et al. 1999). Many of these claims depend on the assumed R_☉ and are normally measured using the data within the solar annulus. In their studies of masers, McMillan & Binney (2010) find that the angular velocity is constrained better ranging between 29.9 and 31.6  km s⁻¹ kpc⁻¹. As a summary of all these works, we assign a prior with a uniform distribution of

The range in ω_LSR of [23, 34] km s⁻¹ kpc⁻¹ corresponds to a range in v_LSR of [196, 289] km s⁻¹ at R_☉ = 8.5 kpc.

We have the σ_r run for two different tracers, i.e., giant and BHB stars, labeled with red and black in Figure 1. While the σ_r run of the giant data is measured out to a larger distance (r ∼ 155 kpc) than of the BHB stars (r ∼ 40 kpc), the giant stars have comparatively larger uncertainties in r and σ_r(r) than BHB stars. Importantly, the β(r ⩽ 25 kpc) profile that goes into our modeling is also unknown for the giant data but known for the BHB data. It is therefore more sensible to take the best portion of the data in hand and combine them. Therefore, for our final round of measurements, we take an adjoined σ_r(r): the BHB data from K12 in the range 12 ≲ r/kpc ≲ 40 and the giant star data in the range 40 ≲ r/kpc ≲ 155. Note that our model for the σ_r(r) (see Figure 3) does not predict the flattening out of the profile in the inner region r ≲ 12 kpc. However, the first two σ_r(r) values in the figure for the BHB sample show a clear flattening. This is something that needs to be investigated in the future; currently we ignore these data points.

During an MCMC run, for every proposed set of values of model parameters, particularly one defining the potentials, we get a prediction for ω_LSR (shown in Figure 5 by a hatched histogram). We are able to constrain the ω_LSR = 30.2 ± 1.2 km s⁻¹ kpc⁻¹. Interestingly, this is within the uncertainty range of the Reid & Brunthaler (2004) result 28.8 ± 0.2 km s⁻¹ kpc⁻¹ for V_☉ = 12.24 km s⁻¹ obtained using the Sgr A*'s proper motion. Also, our estimate falls in the prescribed range 29.9–31.6 km s⁻¹ kpc⁻¹ in McMillan & Binney (2010), who use the proper motions of masers. Our measurement therefore can be taken as an independent measurement of the v_LSR at R_☉ = 8.5 kpc. Since our constraint of ω_LSR has larger uncertainties compare to the one obtained using Sgr A* data, from here on unless otherwise mentioned we use $\mathcal {N}(\omega _{{\rm LSR}}| 28.8, 0.2)$ as a prior instead of a uniform distribution.

Zoom In Zoom Out Reset image size

In this paper we use an Miyamoto–Nagai (MN) disk and a spheroidal bulge. The properties of the disk and the bulge are mainly governed by the terminal velocity data. They are also quite sensitive to the prior on v_LSR, which in turn depends on the chosen value of R_☉. Overall the disk mass is around 10¹¹ M_☉ and the bulge mass is around 10¹⁰ M_☉. The structural parameters like b and s are difficult to measure. The bulge mass has a mild dependence on s, but other than that, b and s have little effect on other parameters (see Figure 6). We find that a+b is well constrained by the data. This can be seen by the strong anticorrelation between them and a narrow spread around it. Our inability to measure b is due to the following two reasons. First, the priors from terminal velocity data and Sgr A* we use basically provide v_circ(R)|_{z = 0}, and this is sensitive only to the sum a+b. Second, in our setup the halo kinematics responds to forces averaged in radial shells, which decreases sensitivity to b.

Since b and s are not well determined by the data, the choice of prior becomes important for them. Traditionally, a double exponential form is used to fit the disk density. By fitting exponential disks to mono-abundance populations, Bovy et al. (2012) finds scale lengths to be roughly in the range 2–4.5 kpc and scale heights to be roughly in the range 0.2–1 kpc. Here we use an MN disk, so to facilitate comparison we derive the appropriate scaling factors. If we fit the surface density of an MN disk with an exponential form, for 2 < R/kpc < 8.5, we get a scale length of about 0.82a. Fitting the density in the vertical direction, in the range 0.5 < z/kpc < 2.0 near the Sun, we get a scale height of 1.75b. We adopt a uniform prior for b in the range 0–0.5 kpc, which is within the range expected for exponential disks. The value of a that we get is also within the range of expectation. For s we adopt a uniform prior in the range 1–3 kpc, which is around the value 1.9 kpc as suggested by Binney & Tremaine (2008).

Our best-fit estimate for the anisotropy in outer parts β_s is 0.4 ± 0.2. It is interesting that we are able to constrain β_s; the reason is as follows. In the innermost region, the terminal velocity data and prior on ω_LSR provide information about the bulge and also to some extent disk parameters. In the region 12 < r/kpc < 25 the BHB anisotropy is already known and the kinematics when put in the Jeans equation provides estimates for r_b, M_vir, c, and disk parameters. Now, beyond r > 25 kpc where β_s is introduced, it is in some sense the only unknown.

We estimate the mass of the dark matter halo to be $M_{{\rm vir}}=0.80^{+0.31}_{-0.16} \times 10^{12} \, M_{\odot }$ and the concentration $c=21.1^{+14.8}_{-8.3}$. Corresponding values for the virial radius R_vir and virial velocity are found to be $239.1^{+27.6}_{-16.6}\ {\rm kpc}$ and $120.2^{+13.9}_{-8.3}\ {\rm km\ s^{-1}}$, respectively. It can be seen in Figure 8 that there is a strong anticorrelation between c and M_vir. The upper bound on c is not as well constrained as the lower bound. Simulated virialized halos in ΛCDM cosmology, in general, predict an inversely proportional relation between their mass and concentration (e.g., Bullock et al. 2001; Macciò et al. 2007; Duffy et al. 2008; King & Mead 2011). The dashed line in Figure 8 shows one such relation, $c = 327.3\, M_{{\rm vir}}^{-0.12}$, adopted from Macciò et al. (2007). We see that the prediction of the ΛCDM simulation for the dark matter halo in the range 10¹¹ ⩽ M_vir/M_☉ ⩽ 10¹³ passes through our measurements. However, note that the predictions are for pure dark matter simulations and do not include baryonic processes such as cooling, star formation, and feedback. The collapse of gas due to cooling leads to adiabatic contraction of the dark matter halo, which increases its concentration. Feedback, on the other hand, can have the reverse effect. Therefore, it is difficult to comment whether the concentration we get is typical or atypical of the Milky-Way-sized galaxies.

Zoom In Zoom Out Reset image size

To answer this, we now run the MCMC separately over a subset of our sample of the halo giant and the BHB samples taken from Figure 1. We select BHB and giant stars in a common radial range, i.e., r ≲ 40 kpc. For the BHB sample we estimate $c = 22.8^{+12.1}_{-7.8}$ and $M_{{\rm vir}}= 0.74^{+0.18}_{-0.12} \times 10^{12} \, M_{\odot }$, whereas for the giant sample we estimate $c = 20.4^{+13.8}_{-8.9}$ and $M_{{\rm vir}}= 0.90^{+0.46}_{-0.26} \times 10^{12} \, M_{\odot }$. The c − M_vir joint-likelihood distributions for the above two data sets are shown by the blue and green contours in Figure 8, respectively. The distributions for both samples seem to be in good agreement, except for the fact that it is more puffed for giants than the BHB sample. We find that this is mainly due to larger uncertainties in the σ_r(r) values for the giants in comparison with the BHB sample, which we verified by swapping the error bars in BHBs and giants. To conclude, the kinematics of the two halo populations, namely, BHBs and giants, are consistent with the fact that they both feel the same Galactic potential, at least within the radius r ≲ 40 kpc.

Kinematics of halo stars allow us to constrain the density profile of the halo. We had modeled the density profile in the inner region by a double power law with fixed slopes, but the location of the break radius was kept free. For BHBs we measure r_b to be $17.1^{+1.2}_{-1.0} \ {\rm kpc}$, and for giants we measure it to be $22.1^{+4.1}_{-3.1} \ {\rm kpc}$. A reason for the different r_b for these two halo populations could be the uncertainties in the distances, which are larger for the giants than for the BHBs. Our estimate for the break radius is slightly smaller than ∼27 kpc as claimed in Deason et al. (2011) or ∼25 kpc as found in Watkins et al. (2009). Interestingly, our estimates are in good agreement with the recent study of RR Lyrae stars in Sesar et al. (2013), who suggest a break in the power law at a much smaller radius of ∼16 kpc. A smaller break radius also complies with the study of SDSS main-sequence turnoff stars in Bell et al. (2008), who conclude that the slope of the density profile at r ≲ 20 kpc should be shallower in comparison with the radius outside this range. Here, a point worth noting is that our estimate of the break radius is linked to the kinematic features, whereas all the above values from the literature are inferred from the studies of the spatial distribution.

To find an answer, we now run the MCMC hammer over the halo giant stars spanning 5 ≲ r/kpc ≲ 155. The red contour in Figure 8 shows the corresponding c − M_vir joint distribution, which is found to almost coincide with the green contours obtained for a giant catalog with r ≲ 40 kpc. Similarly, a comparison of the blue contour in the figure, which is for the BHB sample within 12 ≲ r/kpc ≲ 40, with the black contour, which is for a combined sample of BHB and giant stars in distance ranges 12 ≲ r/kpc ≲ 155, shows that they are similar. This suggests that the giant data in the range 40 < r/kpc < 155 do not add much to our knowledge of c and M_vir, i.e., the potential of the Galaxy. This is mainly due to our adoption of a parameterized form for the distribution of dark mater, namely, the NFW profile. The NFW profile predicts that for r ≫ R_vir/c the density falls as r⁻³. Therefore, data that extend out to about two times the scale radius R_vir/c should be sufficient to constrain the two independent parameters M_vir and c of the NFW profile. However, if one wants to really compute the density out to the virial radius, e.g., by nonparametric schemes, then kinematic data till the virial radius would be required. In our case the distant giant stars (r > 40 kpc) turn out to be useful to probe the density distribution of the tracer population, i.e., the stellar halo. We find $r_t = 97.7^{+15.6}_{-15.8}\ {\rm kpc}$ and $\Delta = 7.1^{+7.8}_{-4.8}\ {\rm kpc}$. It is interesting to note that the hydrodynamical simulations (e.g., Abadi et al. 2006) to investigate the properties of luminous halos surrounding isolated galaxies do not predict a truncated halo; rather, they find that the halo extends to the virial radius. In the future, our results regarding the density profile of the outer stellar halo should be useful for testing theories of stellar halo formation. Recently, Deason et al. (2014), using A-type stars from SDSS, find a sudden drop in density profile of the stellar halo as traced by BHBs and blue stragglers, lending further support to our kinematics detection of such a drop.

The number of subhalos of a given mass scales directly with the host halo mass (Springel et al. 2008). Therefore, an accurate estimate of the Galaxy mass has importance in understanding the missing satellite problem. One interpretation of the problem (Kauffmann et al. 1993; Klypin et al. 1999; Moore et al. 1999; Bullock 2010) is that the ΛCDM paradigm predicts a larger number of massive subhalos for the Milky-Way-sized halo (e.g., Boylan-Kolchin et al. 2011). The problem can be solved if the mass of the Galaxy is low. Figure 5 in Wang et al. (2012) allows us to directly compare the host halo mass with the probability of containing three or less than three subhalos with v_max > 30 km s⁻¹ (maximum value of the circular velocity). The relation is inferred from the studies of the halos obtained from the Millennium Simulation series, Aquarius and Phoenix projects. For a direct comparison we scale our measurement of c, M_vir to compute the mass M₂₀₀ interior to r₂₀₀ from the center of the halo, at which the mean density is 200 times the critical density. We obtain $M_{200} = 0.72^{+0.24}_{-0.13} \times 10^{12} \, M_{\odot }$ and corresponding concentration $c_{200}= 16.2^{+11.6}_{-6.7}$. Figure 5 in Wang et al. (2012) suggests that for the mass equal to our M₂₀₀ there is ∼70% probability that the Galaxy should host three or less than three subhalos with v_max > 30 km s⁻¹. Interestingly, from observations it is known that there are only three brightest satellites of the Galaxy, namely, the Small Magellanic Cloud, Large Magellanic Cloud, and Sagittarius dwarf have v_max ∼ 30 km s⁻¹. This at least suggests that there is no discrepancy between the observed number of luminous satellites with v_max > 30 km s⁻¹ and the number predicted by ΛCDM. Furthermore, Vera-Ciro et al. (2013) also concludes that for Milky Way mass ∼8 × 10¹¹ M_☉, which is similar to our mass, the number and internal dynamics of the dwarf spheroidal satellites of our Galaxy are consistent with the predictions of the ΛCDM model. Therefore, we remark that the scarcity of massive subhalos is not a failure of the ΛCDM paradigm but a repercussion of assuming higher virial mass for the Galaxy.

Another impact of the Galaxy mass is in describing the overall dynamics of the orbiting satellite galaxies. For our low estimate of Galaxy mass, are the satellites still bound is a natural question to ask. To study this, we measure the escape velocity v_esc using

where

Our estimate of v_esc as a function of the galactocentric radius r is shown in Figure 9. The stars in the figure show the total velocities for the named Milky Way satellite galaxies. The velocities are computed from a recent compilation tabulated in Table 1 of Pawlowski & Kroupa (2013). Also, the velocities are corrected for our assumption of the velocity of the local standard of rest, i.e., 245 km s⁻¹ at R_☉ = 8.5 kpc, except Leo II, which seems to be marginally unbound. We can conclude from Figure 9 that all the given satellites are bound despite our low estimate for the Galaxy mass.

Zoom In Zoom Out Reset image size

In Figure 10 we present the cumulative mass M(< r) of the Galaxy. It is computed using the formula for the centrifugal equilibrium

where dΦ/dr is taken from Equation (10) and uses spherical averaging. Infalling satellites are destroyed by their host's gravitational potential, resulting in tidal streams. Attempts to model these streams (Newberg et al. 2010; Law & Majewski 2010; Carlin et al. 2012) also provide an alternate constraint on the Galaxy mass. In Figure 10 the black dot with error bar is the Law et al. (2005) estimate of the Galaxy mass within 50 kpc obtained by modeling the Sagittarius dwarf spheroidal galaxy tidal streams. It is interesting to note that our estimate given the range of uncertainty agrees well with this result. However, the diamond point, which is mass within r < 60 kpc obtained by Newberg et al. (2010) by modeling the Orphan Stream, is significantly smaller than our prediction. One possible reason could be that they model the orbit and not the stream, and the possible misalignment between stream and orbit could bias the result (Sanders & Binney 2013a).

Zoom In Zoom Out Reset image size

The local (at R_☉) dark matter density provides a strong basis for the experimental endeavors for indirect detection of the dark matter; see Strigari (2013) for a review of the topic. Therefore, determination of the local mass distribution, the work originally pioneered by Oort (1932), has recently received a great deal of attention. Our best-fit model of the halo potential allows us to compute the local dark matter density $\rho ^{\rm DM}_{\odot }$, which we measure to be $\rho ^{\rm DM}_{\odot } = 0.0088^{+0.0024}_{-0.0018} \, M_{\odot } \ {\rm pc} ^{-3}$, equivalent to $0.35^{+0.08}_{-0.07}$ GeV cm⁻³. Our result is in good agreement with the recent estimates of 0.3 ± 0.1 GeV cm⁻³ by Bovy & Tremaine (2012), 0.40 ± 0.04 GeV cm⁻³ by McMillan (2011), or 0.389 ± 0.025 GeV cm⁻³ by Catena & Ullio (2010). However, we note slightly lower estimates of 0.007 M_☉ pc⁻³ given in Holmberg & Flynn (2000), who utilize a catalog of A–F stars obtained from the Hipparcos survey and of 0.0065 ± 0.0023 M_☉ pc⁻³ given in Zhang et al. (2013), who utilize K dwarf stars from SDSS/SEGUE.

Also, we measure the local escape velocity, using Equation (15), to be $v_{{\rm esc},R_\odot } = 550.9^{+32.4}_{-22.1}\ {\rm km\ s^{-1}}$. Our estimate seems to be slightly higher but within the range of uncertainties of $544^{+64}_{-46}\ {\rm km\ s^{-1}}$ found using the high-velocity halo stars in Smith et al. (2007). Moreover, the most recent estimate of $v_{{\rm esc},R_\odot }$, again using the high-velocity stars, is provided to be $544^{+64}_{-46}\ {\rm km\ s^{-1}}$ (Piffl et al. 2014). There $v_{{\rm esc},R_\odot }$ is defined to be the minimum speed required to reach three virial radii, where R_vir = 180 kpc. For a fair comparison we redefine Equation (15) to be equal to $\sqrt{2|\Phi (R) - \Phi (3 R_{{\rm vir}})| }$ and compute $v_{{\rm esc},R_\odot } = 528^{+24}_{-17}\ {\rm km\ s^{-1}}$ for R_☉ = 8.5 kpc, which is in the lower range of quoted values in Piffl et al. (2014).

Yet another quantity of interest is whether our disk is maximal (Carignan & Freeman 1985; van Albada et al. 1985). A convention (Sackett 1997; Courteau et al. 2014) is that in a maximal disk 72% of the total rotational support $v_{\rm circ} ^{\rm total}$ is contributed by a disk $v_{\rm circ} ^{\rm disk}$, i.e.,

where R_max is a radius at which $v_{\rm circ} ^{\rm disk}(R)$ is maximum. For our model, we find R_max = 7.4 kpc and F = 0.5 or $v_{\rm circ} ^{\rm disk}(R_{\rm max})/v_{\rm circ} ^{\rm total}(R_{\rm max}) =0.7$, i.e., at this radius 30% of the total rotational support is by a disk.

Because of a lesser contribution of the disk to the total rotation curve, we find the Galaxy disk at R_max to be submaximal. Recently, Bovy & Rix (2013), using SEGUE dwarfs, found a slightly higher value of F = 0.69 and concluded the disk to be maximal. They assumed an exponential disk, for which F is measured at R = 2.2 R_d (scale length of the exponential disk). It is, however, interesting to note that in their studies of 12 out of 15 distant spiral galaxies, Kregel et al. (2005) find on average F = 0.28, which is nearly half of our estimate for the Galaxy.

Finally, we check whether the local (at R_☉) column/surface density of our best-fit model is within an expected range. The surface density is computed using

The contribution of the disk to the surface density is strongly dependent on our prior for disk scale length b, so it is not really a prediction of our analysis. We find that $\Sigma ^{\rm total}_{\odot }$ is $94.0^{+16.6}_{-20.3} \, M_{\odot } \ {\rm pc} ^{-2}$. This is slightly higher than 71 ± 6 M_☉ pc⁻² obtained by Kuijken & Gilmore (1991) and 74 ± 6 M_☉ pc⁻² obtained by Holmberg & Flynn (2004), but within the range of large uncertainties we have. We note that recent measurements by Bovy & Rix (2013), Zhang et al. (2013), who use the kinematics of the SDSS/SEGUE dwarfs, suggest a slightly lower value of $\Sigma ^{\rm total}_{\odot } = 68 \pm 4 \, M_{\odot } \ {\rm pc} ^{-2}$.

The results presented here are potentially subject to systematic uncertainties such that the reader should be cautioned. We use fiducial isochrones for the distance estimation of the giants. But there might be systematics associated with them, and this will have an effect on our distance estimates. An independent way to validate our distance measurement would be to use an estimator shown in Schönrich et al. (2012), but this requires proper motions. It remains to be seen if any available proper motions of the distant giants are accurate enough for this method to be successfully applicable. Additionally, while measuring the distance to the halo K giants, we do not take into account uncertainties in the reddening estimate. We simply use colors dereddened according to Schlegel et al. (1998) extinction maps. This could introduce systematics in our distance estimate. However, we believe that such systematics, if any, would be insignificant. Firstly, the Schlegel et al. (1998) maps are accurate for high-latitude stars, and 98% of our sample has |b| > 20°. Secondly, although, Schlegel et al. (1998) maps provide total extinction, they are appropriate for halo stars as most of the dust is confined to the disk.

We use a sample of giants from SEGUE, which are subjected to a proper-motion restriction of 11 mas yr⁻¹ (Yanny et al. 2009). This can potentially introduce a systematic bias for nearby giants, but the giants that we use for our main analysis are beyond 40 kpc, and for them the above-mentioned proper-motion limit safely includes all bound halo stars.

We assume that each tracer population is in Jeans equilibrium. In any case, it is important to note that we can determine the Galactic potential only to the extent that the phase-space distribution of tracer stars is in equilibrium (Binney 2013). However, if the tracer population under study is a superposition of multiple populations, then the Jeans equation should be applied separately for each population. In earlier studies with BHB stars (Kafle et al. 2013; Hattori et al. 2013), correlation between metallicity and kinematics of halo targets was found. A similar correlation could exist for giants. So ideally, if one has a larger sample of stars, one should treat metal-rich and metal-poor populations separately. Moreover, close to the disk the potential is not spherically symmetric, so strictly speaking the spherical Jeans equation is not valid. This can potentially bias the estimate mass. This should be investigated in the future.

We model the disk by the MN form. In reality, the disk is described better by double-exponential functions. Moreover, the disk of the Milky Way is probably a superposition of multiple populations, each with different scale height and length. These facts can potentially bias our results. Finally, throughout our analysis we assume R_☉ = 8.5, and systematics due to our adopted value can be expected. Claims for a wide range in R_☉ between ∼7.7 and 8.8 kpc exist in the literature (e.g., Reid 1993; Nikiforov 2004; Gillessen et al. 2013; Francis & Anderson 2014; Reid et al. 2014, etc.). To study the influence of our adopted value of R_☉, we reanalyze our data with R_☉ = 8 kpc. Firstly, we measure the σ_r(r) run, which we find to be similar to one shown in Figure 1. Secondly, we run the MCMC model fitting to the above σ_r(r) run. Assuming R_☉ = 8 kpc means that the prior in v_LSR changes, and this can have a significant influence in our measurement. It is mainly because we do not have halo data within r ≲ 12 kpc, and the information about disk and bulge properties mostly comes from the assumption about v_LSR at R_☉. We find that for R_☉ = 8 kpc our model parameters are M_vir = 1.2 ± 0.3 × 10¹² M_☉, $c = 15.2^{+5.6}_{-3.2}$, $M_{{\rm disk}}= 0.71^{+0.11}_{-0.13} \times 10^{11} \, M_{\odot }$, a = 3.0 ± 0.6 kpc, $M_{{\rm bulge}}= 0.70^{+0.31}_{-.32} \times 10^{10} \, M_{\odot }$, and s = 2.1 ± 0.6 kpc, whereas r_b, r_t, Δ, and β_s remain unchanged.

Recall that our final results shown in Table 2 or Figure 6 are obtained for a combined sample of two halo populations, i.e., BHBs within a radius r ≲ 40 kpc and giants outside r ≳ 40 kpc. As discussed earlier in Section 4.1, the density profile for same sample of BHB stars as ours is already computed in Deason et al. (2011), and this is what we assume in our analysis. In the region where giant star data is used we assume the density profile to be a power law with slope of −4.5, which is consistent with the findings for BHB (Deason et al. 2011) and RR Lyrae stars (Watkins et al. 2009). However, there is no direct measurement of the density of halo giants. Thus, to test the sensitivity of our results to the density profile, we run the MCMC simulation for two different values of slope indices, −4 and −5. From our fits we find that M_vir and c are directly proportional to the density slope. For each linear step of −0.5 from −4 to −5 we find that M_vir increases by ∼30% for each step, whereas c increases by ∼20%. For the same steps from −4 to −5, we find that M_disk and M_bulge decrease by ∼22% and ∼20%, respectively. Quantities such as r_b, r_t, Δ, and disk and bulge scale parameters remain unchanged. Moreover, the shape of the posterior distributions of the parameters and hence uncertainties remain the same for all the cases. It should be noted that the value of β_s remains nearly the same for both the cases, i.e., 0.5 ± 0.2 for −4 and 0.4 ± 0.2 for −5. This is not unexpected, as data in 12 < r/kpc < 25 for which observed β is available are enough to constrain M_vir and c. So any change in density profile affects the enclosed mass. When data beyond 25 kpc are added, the increase in power-law index is compensated by the increase in mass so as to leave β_s unchanged.

Before concluding, it is worth mentioning that we noted the parallel work by Bhattacharjee et al. (2014). It has some similarities with our works, e.g., uses of the Jeans equation, SDSS/SEGUE BHB and giant star catalog, etc. While estimating the mass using the Jeans equation, we note that they use the density profile of the spectroscopic sample, which has a selection bias. In reality, the Jeans equation requires the underlying density profile of the tracers and not of the sample. Lastly, it has been found in the simulations that halo stars, satellites, and the dark matter halo have different orbital properties (Abadi et al. 2006; Sales et al. 2007). Hence, assuming a constant anisotropy for both field stars and satellites could introduce additional systematic uncertainties in their mass estimate. However, it should be noted that they do not model the disk, bulge, and halo separately, but only provide an estimate for the total mass M(< r) within some radius.