KINEMATICS OF THE STELLAR HALO AND THE MASS DISTRIBUTION OF THE MILKY WAY USING BLUE HORIZONTAL BRANCH STARS

Understanding the formation of stellar halos provides vital insight into the formation history and the evolution of galaxies (Majewski 1993; Freeman & Bland-Hawthorn 2002; Helmi 2008). Under the currently favored ΛCDM model of galaxy formation, the stellar halos are thought to have been built up, at least in part, by accretion of satellite galaxies (White & Rees 1978; Searle & Zinn 1978; Helmi & White 1999). Discoveries of structures like the Sagittarius dwarf galaxy (Ibata et al. 1994, 1995; Majewski et al. 2003; Belokurov et al. 2006), the Virgo overdensity (Jurić et al. 2008), the Triangulum–Andromeda structure (Rocha-Pinto et al. 2004; Majewski et al. 2004; Martin et al. 2007), and the low-latitude Monoceros ring (Newberg et al. 2002) lend further support to the idea of the stellar halo being formed by accretion. Other than accretion, a part of the stellar halo could also be formed by in situ stars. Recent hydrodynamical simulations (Abadi et al. 2006; Zolotov et al. 2009; Font et al. 2011; McCarthy et al. 2012) of galaxy formation suggest that in the inner regions, the stellar halo might be dominated by in situ stars, whose kinematic properties are distinct from the accreted stars.

Observational evidence of multi-component halos has been found in dynamical studies of Sloan Digital Sky Survey (SDSS) calibration stars (Carollo et al. 2007, 2010; Beers et al. 2012), in rotational behavior in metallicity bins (Deason et al. 2011a), in kinematics of the Galactic anti-center and North Galactic Pole population (Kinman et al. 2007), in chemical properties (de Jong et al. 2010), and also in the age difference between in situ and accreted halos (Kalirai 2012). There also exists a counterclaim by Schönrich et al. (2011), who demonstrate that the evidence of a retrograde signal in the outer halo in the Carollo et al. (2010) sample is weak and is due to a manifestation of incorrect distance estimates. Investigating whether any signal of a multi-component halo could also be seen in the dispersions of the velocity of the halo population is important. Ultimately, studying the velocity dispersion profiles of the halos and comparing them with simulations might help to constrain the model of galaxy formation.

The lack of proper motions and the slightly off-centered position of the Sun with respect to the galactic center poses a unique challenge in studying the kinematics of the stellar halo. At distances much larger than R_☉, line-of-sight velocity is the same as radial velocity with respect to the galactic center. Hence, the σ_r profile is easy to compute at large distances, a fact that has been well studied; Battaglia et al. (2005) studied the line-of-sight velocity dispersion (σ_los) of a mixture of 240 halo objects and found that σ_los decreases monotonically beyond r ∼ 30 kpc. In the outermost halo at r ∼ 100 kpc, Battaglia et al. (2005) and recently Deason et al. (2012b) both see a significant drop in the σ_los value to ∼50 km s⁻¹.

Conversely, De Propris et al. (2010) studied 666 blue horizontal branch (BHB) stars from the 2QZ Redshift Survey and found that the velocity dispersion profile increases at large distances. However, using 910 distant halo A-type stars, Brown et al. (2010) found that there is a mean decline of -0.38 ± 0.12 km s⁻¹ kpc⁻¹ in σ_r over 15 <r/kpc < 75. More recently, Xue et al. (2008) used 2401 BHB halo stars within 60 kpc and measured a slower decline in σ_los compared to earlier studies. At small r, it has been difficult to derive the σ_r profile, and the only attempt to measure σ_r in the inner halo was undertaken by Sommer-Larsen et al. (1997). They found a sharp fall in σ_r from ∼140 to ∼100 km s⁻¹ at r ≈ 12 kpc, although they assumed the circular velocity to be constant.

In a solar neighborhood, one can get useful proper motions, with which all three velocity dispersions (σ_r, σ_θ, σ_ϕ) can be measured. Smith et al. (2009b) used the full phase space information of ∼1700 halo subdwarfs from the solar neighborhood (<5 kpc) and determined the velocity dispersions to be (σ_r, σ_θ, σ_ϕ) = (143 ± 2, 82 ± 2, 77 ± 2) km s⁻¹. Also, Bond et al. (2010) analyzed the proper motions of a large sample of main-sequence stars within the solar neighborhood (<10 kpc) and found σ_r = 141 km s⁻¹, σ_θ = 75 km s⁻¹, and σ_ϕ = 85 km s⁻¹. A summary of the estimated values of the velocity dispersions found in the literature is given in Table 1.

Due to the lack of reliable proper motions of the halo field stars, even fundamental quantities like the tangential components of the velocity dispersions as well as the anisotropy (β) are still poorly understood beyond the solar neighborhood. The situation is not, however, as hopeless as it might seem. At small r, using only the line-of-sight velocity it is possible to put some constraints on these quantities using maximum likelihood techniques, where a model or a distribution function (DF) needs to be specified a priori. In an analysis of 1170 BHB stars ranging from 5 ≲ d/kpc ≲ 96, where d is now the stellar distance rather than a radius, Sirko et al. (2004b) fit an ellipsoidal distribution of velocities and find that the halo is isotropic. Similarly, Deason et al. (2011a, hereafter D11) fit a constant anisotropy model (power-law distribution function) to 3549 BHB stars constructed from the SDSS Data Release 7 (DR7) and find that the halo between r = 10 and 25 kpc is tangential, whereas the distant halo within 25 < r/kpc < 50 is radial. D11's claims for the tangential inner halo contrast with the result by Bond et al. (2010), who found the inner halo to be radial in the same region (d ≈ 10 kpc). D11 assume the DF to be such that the tracer density and the potential are both power laws. Without a priori knowledge of the density slope, their estimates will have some systematics. Additionally, the potential was also kept constant in their analysis and thus can bias the results due to the degeneracy between the potential and the anisotropy.

In their most recent work, Deason et al. (2012a, hereafter D12) allow both the potential and β to vary, and thereby break the degeneracy, and find that the outer halo within 16 < r/kpc < 48 is radial with β = 0.5^+0.1_−0.2. Previous estimates of the velocity anisotropy (β) in the solar neighborhood and the nearby halo are summarized in Table 1. All the above estimates of an anisotropy of the distant halo are done in a large radial bin. Their results might be accurate for the given radial bins and could also be the actual anisotropy of the halo, provided the anisotropy remains nearly constant or monotonic throughout. On the other hand, if the actual β of the system is not monotonic but has a more complex shape, then estimating it in the larger bin will just capture an overall property of that bin.

Theoretically, there are families of the DF of the Henon (1973) type that result in a constant anisotropic system, as well as families of models that have their own anisotropy profiles, which include the Osipkov–Merritt model (Osipkov 1979; Merritt 1985a, 1985b), Gerhard (1991), Cuddeford (1991), Baes & van Hese (2007), and few more with the Hernquist potential-density model in Baes & Dejonghe (2002). The question we ask is how well do these anisotropy profiles match the anisotropy of the Galaxy? More fundamentally, how well do we know the anisotropy of the halo? To this end, we compute the beta profile with much finer spatial resolution and without any prior assumptions about the density or the potential.

Another use of studying the kinematics of the stellar halo is to constrain the mass and the potential of the Milky Way. Most of the methods to estimate the mass require knowledge of the anisotropy parameter β. Without the unbiased estimate of the velocity anisotropy, constraining the mass of the system via the Jeans equation could be uncertain by 73% (for −4.5 < β < 0.44), as found by Watkins et al. (2010) in their estimates of the mass of the Galaxy. Several other authors have also used this assumption to estimate the mass of the galaxy (Dehnen et al. 2006; Gnedin et al. 2010; Deason et al. 2012b). Precautions need to be taken when making an arbitrary assumption about the anisotropy. It has been found 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 (Battaglia et al. 2005; Dehnen et al. 2006) could introduce systematic uncertainties in the mass estimate. Ideally, in order to break the degeneracy, we must have a separate estimation of the radial velocity dispersion, the velocity anisotropy, and the underlying density of the population, as pointed out by Dehnen et al. (2006).

The orbital evolution of the Magellanic clouds (Lin & Lynden-Bell 1982; Besla et al. 2007), the local escape speed (Smith et al. 2007), the timing argument (Li & White 2008), and the study of the kinematics of the tracers population (Kochanek 1996; Xue et al. 2008; Gnedin et al. 2010; Watkins et al. 2010) are the methods undertaken in order to constrain the mass of the Galaxy. Summarizing all these attempts to constrain the mass of the Galaxy, the consensus can be found between 0.5–3.5 × 10¹² M_☉. Recently, using BHB stars, Xue et al. (2008) estimate the mass of the Milky Way to be 0.91^+0.27_−0.18 × 10¹² M_☉. However, they make an assumption that the variation of (v_circ/v_los), the ratio of the circular to the line-of-sight velocity, with radius is the same as in simulations. In this paper, we estimate v_circ as far out as possible without any assumption and then use it to estimate the dark matter halo mass.

This work focuses mainly on the study of the kinematics of the stellar halo in order to present an unbiased estimation of the velocity dispersions, anisotropy parameter, and circular velocity as a function of radius to the extent the data will support. We use a DF which does not require any assumptions to be made a priori about the density profile or the potential. We then use our measurements of velocity dispersions to estimate the rotation curve of the Galaxy. With the disk and bulge mass already constrained from Sofue et al. (2009), we focus on constraining the dark matter halo mass. Using the circular velocity curve (v_circ(r)), we can estimate β(r) out to as far as σ_r(r) is available. Finally, we compare our results with simulations in which the halo is formed purely by accretion.

This paper is organized as follows. In Section 2, we discuss the theoretical aspect of our analysis, the methodology adopted, and the details about the sample. In Section 3, we present our results for the velocity dispersions, compare the results between the alternative models, and investigate the contribution of the halo substructures. In Section 5, we present our estimation of the mass of the Galaxy. Results are then compared with the simulations in Section 6. In Section 7, we present our conclusion and discuss the results.

We are interested in calculating the velocity dispersions (σ_r, σ_θ, σ_ϕ) for the stellar halo. However, the data we have are line-of-sight velocities. To proceed, we need to make some assumptions about the position and the velocity of the Sun with respect to the galactic center. We assume R_☉ = −8.5 kpc, the velocity of the local standard of rest (LSR), v_LSR, is taken to be an IAU adopted value = 220 km s⁻¹, and the solar motion with respect to LSR (U, V, W)_☉ = +11.1, +12.24, +7.25 in km s⁻¹ (Schönrich et al. 2010). Spherical and heliocentric coordinate systems are expressed in terms of (r, θ, ϕ) and (d, l, b), respectively.

2.1. Distribution Function

The DF, f, is defined such that f(x, v) d³x d³v is the probability of finding a randomly picked star in a phase-space volume d³x d³v. In general, we consider the stellar halo to be an anisotropic spherical system. The anisotropy is defined as

$$\begin{equation} \beta = 1 - \frac{\sigma _{\theta }^{2} + \sigma _{\phi }^{2}}{2\sigma _{r}^{2}}, \end{equation} \tag{ 1 }$$

where σ_r, σ_θ, and σ_ϕ are the velocity dispersions in spherical coordinates, and it describes the orbital structure of the system. The values of this parameter range from –∞ for purely circular trajectories to 1 for purely radial orbits.

Families of DFs that generate the collisionless anisotropic spherical systems with constant or varying velocity anisotropy can be found in detail in Binney & Tremaine (2008). One such DF with a constant anisotropy is

$$\begin{equation} f(E,L) = f(E) L^{-2\beta }. \end{equation} \tag{ 2 }$$

Here, E = Φ(r) − (v²/2) is the relative energy per unit mass and L is the modulus of the angular momentum vector per unit mass. Recently, a DF given by Equation (2) with the energy term from Evans et al. (1997),

$$\begin{equation} f(E,L)\propto E^{(\beta (\gamma -2)/\gamma)+(\alpha /\gamma)-3/2} L^{-2\beta }, \end{equation} \tag{ 3 }$$

was used by Deason et al. (2011a, 2012a) to study the rotation, anisotropy, and mass of the Galactic halo. The parameters α and γ are the logarithmic slopes of the density (ρ∝r^−α) and potential (Φ∝r^−γ), respectively. Hereafter, we refer to this function as D11 DF.

If one is interested in deriving the dispersion profiles, then a simple DF that one can use is the Gaussian velocity ellipsoidal DF (GVE DF). The GVE DF has been used in the context of the stellar halo by Frenk & White (1980) using globular clusters as tracers, and by Sirko et al. (2004b) and Smith et al. (2009a) using halo stars. A GVE DF with rotation about the z-axis is given by

$$\begin{equation} f(r, {\bf v}) = \frac{\rho (r)}{(2\pi)^{3/2}\sigma _{r}\sigma _{\theta }\sigma _{\phi }} {\rm exp}\left[\! {-}\frac{1}{2} \!\left(\! \frac{v_{r}^{2}}{\sigma _{r}^{2}}+ \frac{v_{\theta }^{2}}{\sigma _{\theta }^{2}}+\frac{(v_{\phi }-v_{\rm rot})^{2}}{\sigma _{\phi }^{2}} \!\right) \!\right]. \vspace*{4pt} \end{equation} \tag{ 4 }$$

The DF as given by Equation (4) assumes that the velocity ellipsoid is perfectly aligned with the spherical coordinates, but in general the velocity ellipsoid can have a tilt. Using halo subdwarf stars, Smith et al. (2009b) and Bond et al. (2010) have found that the tilt is small and consistent with zero. Hence, it is safe to ignore the tilt while computing velocity dispersions.

2.2. Parameter Estimation

The information on the proper-motion of the stars in the stellar halo beyond the solar neighborhood (r ≳ 10 kpc) is not accurate enough to properly constrain the tangential motions. Nevertheless, our position in the Galaxy still makes it possible to constrain these quantities by utilizing the tangential information carried by the line-of-sight velocities of the stars. However, for that, we need to marginalize the DF over the unknown quantities, which in this case are the tangential components (v_l, v_b). The expression for the marginalized DF is given by

$$\begin{equation} F(l, b, d, v_{\rm los}|\sigma _{r},\sigma _{\theta },\sigma _{\phi },v_{\rm rot}) = \int\int f {d}v_{l} {d}v_{b}. \end{equation} \tag{ 5 }$$

We use the maximum likelihood method to estimate the model parameters. The log-likelihood function which we maximize is given by

$$\begin{equation} \mathcal {L} (l, b, d,v_{\rm los}|\sigma _{r},\sigma _{\theta },\sigma _{\phi },v_{\rm rot}) = \sum _{i}^{n}{\rm log}F(l_{i}, b_{i}, d_{i}, v_{{\rm los}_{i}}), \end{equation} \tag{ 6 }$$

where n is the number of stars in the system under study. We use Markov Chain Monte Carlo (MCMC) with the Metropolis–Hasting algorithm to obtain the posterior distribution. We quote the central values of the velocity dispersions (σ_r, σ_θ, σ_ϕ) as our initial estimates with 16 and 84% error. Note that for the GVE DF, the density term ρ(r) in Equation (4), is not a function of model parameters and hence does not have any effect on the likelihood analysis.

Once the radial velocity dispersion σ_r and the anisotropy parameter β are evaluated, the Jeans equation (Jeans 1915) can be used to estimate the circular velocity v_circ of the spherical system in equilibrium using the relation

$$\begin{equation} v_{\rm circ}^{2}(r) = -\sigma _{r}^{2} \left[\frac{{d\ln }\rho }{{d\ln } r}+\frac{{d\ln }\sigma _{r}^{2}}{{d\ln }r}+2\beta \right], \end{equation} \tag{ 7 }$$

where ρ∝r^−α is the density of the tracer population, which implies dln ρ/dln r = −α. Throughout the analysis, we assume the density to be a double power law with α = 2.4 (r ⩽ 27 kpc) and α = 4.5 (r > 27 kpc), in agreement with the recent works by Deason et al. (2011b) and Watkins et al. (2009).

For systems with a constant anisotropy and a given v_circ, the solution to the differential equation (7) subject to the boundary condition lim_{r → ∞} σ²_r = 0 reads

$$\begin{equation} \sigma _{r}^{2}(r) = \frac{1}{r^{2\beta }\rho (r)} \int _{r}^{\infty } dr^{\prime } r^{\prime 2\beta } \rho (r^{\prime }) ({d}\Phi /{d}r^{\prime }). \end{equation} \tag{ 8 }$$

Assuming that density and anisotropy are known, we can use this solution to estimate σ_r as a function of r.

2.3. DATA: BHB Stars

Since they are luminous and have nearly constant magnitude, BHB stars are ideal for studying the stellar halo, and this is what we use in our study. We use the BHB catalog published by Xue et al. (2011, hereafter X11) for our analysis. The catalog is composed of 4985 BHB stars obtained from SDSS DR 8 (Aihara et al. 2011). The stars were selected by imposing limits on color and Balmer line profile measurements. Imposing limits on Balmer line profile measurements allows one to remove the main-sequence stars and Blue Stragglers. Further details on BHB candidate selection can be found in Xue et al. (2008) and references therein. To avoid contamination from the disk stars, we restrict our analysis to stars with a distance of |z| > 4 kpc from the galactic midplane. As mentioned earlier, due to the way our likelihood function (Equation (5)) is laid out, this cut in distance above the plane will not introduce any bias. No velocity limits have been imposed to obtain the sample and, thus, for the purpose of kinematic studies, the population of BHB stars we select can be considered to be kinematically unbiased.

For the stars that we study, the angular position is known very accurately, but the distance and radial velocity have some uncertainty associated with them. To get more accurate distances, we recalibrate the X11 distances using a color–magnitude relation derived for the same population from Deason et al. (2011b). The estimated dispersion in g-band magnitudes is 0.13, equivalent to a distance uncertainty of 6%. For the SEGUE (Yanny et al. 2009) radial velocity measurements, 94% of our sample have an uncertainty of less than 8 km s⁻¹.

The galactocentric radial distribution of the final BHB samples is shown in Figure 1. One can see that the distribution peaks at around 16 kpc. Most of the stars are found to lie in the range 10 < r/kpc < 25.

Zoom In Zoom Out Reset image size

We study the kinematics of the halo in radial bins to obtain the radial profile of the model parameters (σ_r, σ_θ, σ_ϕ) and also β. Using only the line-of-sight velocity information, the tangential components, σ_θ and σ_ϕ, are difficult to constrain except in the very inner regions of the halo. However, σ_r can be well constrained in both the inner and outer halos. This means that a relatively larger number of stars (>1000) per bin are required to estimate σ_θ and σ_ϕ as compared to σ_r. Given that we only have about 4000 stars, this means that we cannot measure the σ_θ and σ_ϕ profiles with sufficient spatial resolution. Hence, we employ two different binning schemes or estimators, one for radial velocity and the other for tangential velocity. The estimators are the equi-populated estimator (EPE) and the central moving estimator (CME). In EPE, the data are binned radially with each bin containing an equal number of particles, and this is used for computing σ_r(r). In CME, a set of equi-spaced positions in r are chosen and then at each position an equal number of points on either side of the chosen central value are used to estimate the desired quantity. We use the CME for computing σ_θ(r) and σ_ϕ(r). The crucial difference between the two schemes is that while the bins are non-overlapping in the former, they can be overlapping in the latter. In EPE, the spacing between the bins is directly proportional to the number of particles in each bin n_bin. Hence, if the desired quantity can be estimated with sufficient accuracy by employing small n_bin, then EPE is the desired method. However, if this is not the case, then it is better to use the CME method as the spatial resolution is not directly dependent on n_bin.¹

Finally, for our data, the number density of points in r is highly non-uniform, and hence it is not accurate to assume that the desired quantity has been estimated at the center of the bin. To alleviate this number density bias, for both schemes we compute the final position of the bin as the mean r of the points in the bin.

3.1. Radial Velocity Dispersion Profile (σ_r(r))

Here, we focus on the nature of the σ_r(r) profile of the halo for which we adopt the EPE method with n_bin = 400. As explained previously, σ_r can be measured out to the extent of the data (r ∼ 60 kpc). The values of σ_r obtained from the likelihood analysis are given in Figure 2 and the error bars represent the 1σ confidence interval and are determined from the likelihood fitting.

Zoom In Zoom Out Reset image size

We find that the radial velocity dispersion, σ_r, at the Sun's position (R_☉ = 8.5) is 145.6 km s⁻¹. However, beyond the solar neighborhood, σ_r sharply decreases until r ∼ 15 kpc, after which it decreases much slowly and approaches a value of ∼100 km s⁻¹ at 56 kpc. The error bar in σ_r for the inner halo population is large and is mainly due to the fact that v_los contains less radial velocity information compared to outer parts. The overlaid black dashed line is the linear approximation (≈111–0.31r) for the σ_los profile from Xue et al. (2008). Additionally, other previous attempts to fit the profile for σ_r in the outer parts of the halo (Battaglia et al. 2005; Gnedin et al. 2010; Brown et al. 2010) also found profiles similar to Xue et al. (2008), with a slightly varying slope and normalization. All of these profiles are reasonable estimates of σ_r for the outer halo (d ≫ R_☉), where the assumption σ_r ≈ σ_los holds. In the inner halo (r ≲ 15 kpc), however, the approximation breaks down and σ_r strongly deviates from σ_los. One can see that the deviation of σ_r from σ_los increases as one approaches the center, and at R_☉ the deviation is as high as ∼40 km s⁻¹.

Sommer-Larsen et al. (1997) provide a functional form for fitting the σ_r profiles, which is given by

$$\begin{equation} \sigma _{r}^{2} = \sigma _{0}^{2} + \frac{\sigma _{+}^2}{\pi }\left[\frac{\pi }{2}-{\rm tan}^{-1} \left(\frac{r-r_{0}}{l}\right)\right]. \end{equation} \tag{ 9 }$$

This has a shape that is similar to our σ_r(r) profile and we fit it to find σ₀ = 94.5 km s⁻¹, σ₊ = 122.3 km s⁻¹, r₀ = 13.2 kpc, and l = 2.6 kpc. In this function, the fit parameter σ₀ gives the asymptotic value that σ_r achieves in the outer halo, whereas (σ²₀ + σ₊²)^1/2 gives the approximate value for σ_r in the inner halo. The fit parameters r₀ and l determine the turnoff point and the steepness of the transition of the profile, respectively. A lower value of l gives a steeper transition; this can be seen from the comparison between the Sommer-Larsen fit (l = 7.5 kpc) and the red line in Figure 2, which is our fit with a smaller l.

3.2. Tangential Velocity Dispersion and Anisotropy Profiles (σ_θ(r), σ_ϕ(r), and β(r))

There have been few attempts to constrain the tangential velocity dispersions and most of the studies are either restricted to the solar neighborhood (Smith et al. 2009a; Bond et al. 2010) or for the overall system (Sirko et al. 2004a). The tangential velocity dispersions not only provide information about the anisotropy of the stellar population (through Equation (1)), but together with σ_r also help to measure the mass distribution of the Galaxy.

We estimate the tangential velocity dispersions (σ_θ and σ_ϕ) using CME with n_bin = 1200 stars. Given the quality of the data, we are only able to measure σ_θ and σ_ϕ out to ∼25 kpc. Our estimates of σ_θ and σ_ϕ are shown in Figures 3(b) and (c) by the black dots with error bars. For uniformity, we also estimate σ_r with this binning scheme, which is shown (only out to ∼25 kpc) in Figure 3(a) by the black dots with error bars.

Zoom In Zoom Out Reset image size

In general, the tangential components σ_θ and σ_ϕ near the solar neighborhood are found to be comparatively lower (σ_θ = 85⁺⁸₋₉ km s⁻¹, σ_ϕ = 95⁺⁸₋₈ km s⁻¹) than the radial dispersion σ_r. Figures 3(b) and (c) show that there is a sharp rise in the values of σ_θ and σ_ϕ at r = 17 kpc. Beyond this, σ_ϕ falls, whereas σ_θ rises; given the large uncertainties and low number of independent bins, it is unclear if this is a real or a spurious trend.

By substituting the estimates of the tangential and the radial velocity dispersions obtained using CME (n_bin = 1200) from the above analysis into Equation (1), we compute the corresponding values of the anisotropy constant in the respective bins. As shown in Figure 3(d), the halo within 12 kpc has β ∼ 0.5, whereas the outer halo beyond the turnoff point is nearly isotropic. We discover a significant drop in the β profile at r = 17 kpc. Here, the halo is strongly tangential with β = −1.2. We later confirm that the trend observed in anisotropy is neither due to the manifestation of the systematics introduced by the chosen model (Section 3.3) nor due to the presence of the halo substructures (Section 3.4). It is also shown in Appendix B that the assumed v_LSR and R_☉ have negligible effects upon these estimates. The probable reasons for this sudden turnover in the properties of the halo are discussed in the conclusion.

We know that the consecutive CME bins overlap in radius, and thus the dispersion profiles demonstrated in Figure 3 are a smoothed version of the actual dispersion profiles of the halo. However, to check for any systematic associated with the choice of the binning scheme, we also estimate σ_r, σ_θ, σ_ϕ, and β in traditional equi-populated (EPE) bins (n_bin = 700). The measured values in these bins are shown by the diamond points in Figure 3. If the number of stars per bin is fewer than 700, then it is difficult to constrain σ_θ and σ_ϕ. Even with n_bin = 700, we were only able to constrain the tangential motion until 16 kpc (first three diamond points). Hence, we construct the last bin by grouping all the stars beyond 16 kpc into one bin (rightmost diamond point). More importantly, with this binning scheme, we are only interested to see whether we obtain the corresponding dip or rise (depending on the parameter of interest) seen in Figure 3 or not. We find that except for the rightmost diamond points, all of the other diamond points in the figure are in agreement with our previous estimates of the dispersions (given by black dots). The rightmost diamond points are calculated in a huge bin with more than 50% of the total sample. Particularly for σ_θ, σ_ϕ, and β, given their non-monotonic trend, we do not claim that the last diamond point is the correct estimate of anisotropy at r ∼ 35 kpc.

None of the uncertainties quoted in the above estimates of σ_r, σ_θ, and σ_ϕ include the uncertainties in distances and radial velocities. As mentioned earlier, the errors in distance and radial velocity are quite small, and convolving the model (Equation (4)) with the error functions should not change the results.

3.3. β(r) from Fitting Distribution Function (f(E,L))

It would be interesting to see whether the β profile presented above, in particular the dip seen at r = 17 kpc, is an artifact of our chosen model (GVE) or a real inherent feature of the Galactic halo. In order to pursue it, here we explore the effect of the chosen model on the determination of the velocity anisotropy (β). For the comparative study, the alternative model we choose is the D11 DF (Equation (3)). We consider anisotropy (β) and potential (= Φ₀r^−γ) as free model parameters and constrain them using the maximum likelihood method.

First, we compare the theoretical properties of the DFs at our disposal, namely GVE and D11 DF. In Figure 4, we plot the theoretical line-of-sight velocity distributions (LOSVDs) of models along two separate lines of sight. We also show the LOSVDs in Figure 4 for two different distances representing the inner halo (d = 15 kpc, given by red lines) and the outer halo (d = 50 kpc, given by blue lines). For the inner halo, we assume α = 2.4 and assign radially biased anisotropy (β = 0.4), whereas for the outer halo we assume α = 4.5 and assign tangentially biased anisotropy (β = −2.5), in accordance with the Watkins et al. (2009) and Deason et al. (2011b) estimates for α. Note that the density normalization at the break radius (r = 27 kpc) is assumed to be equal. For the assumed constant potential, the solid lines in Figure 4 (both left and right panels) are the LOSVDs obtained by adopting the D11 model and the dashed lines are the LOSVDs of our GVE model. Recall that our model does not take β directly but demands the information of the velocity dispersion components (σ_r, σ_θ, σ_ϕ) individually. To make the LOSVDs obtained from both the models comparable, we estimate σ_r, σ_θ, and σ_ϕ from the set of values of β, α, Φ₀, and γ chosen to obtain LOSVDs of D11 DF. For an assumed potential power law, we use Equation (8) to first calculate σ_r, to put in our model. For an assumed β, substituting this σ_r in Equation (1) gives the corresponding value for σ_t (= $\ssty\sqrt{\sigma _{\theta }^{2} + \sigma _{\phi }^{2}}$). The figure shows that for all four cases, LOSVDs obtained from both models match well. Naively, from these perfect matches of the LOSVDs at different lines of sight and distances, it can be anticipated that the estimation of β(r) with both models should also match.

Zoom In Zoom Out Reset image size

Now, we estimate the β(r) profile using exactly the same sample of BHB stars in the same radial bins as in Section 3.2 (CME, n_bin = 1200) but with D11 DF. Blue points in Figure 5 demonstrate the β profile obtained by fitting D11 DF. Here, we use brute-force grid-based analysis to constrain the model parameters β, Φ₀, and γ. To facilitate the comparison, our estimates from Figure 3(d) are overplotted in Figure 5 and are shown by the black points. Figure 5 shows that within the range of uncertainties, the measured values of β agree with both models (D11 DF and GVE models). However, a slight bias exists, in the sense that β obtained from GVE DF is higher generally. The reason for this discrepancy lies in the fact that β obtained from D11 DF has a dependence on α. Hence, unless the underlying α value of the sample is known exactly, a mismatch is expected. The estimated value of β increases with the adopted value of α (see Figure 3 of D11). This suggests that the actual value of α is even higher than the one adopted here (2.4).

Zoom In Zoom Out Reset image size

In order to estimate the quality of the constraints obtained from the likelihood analysis, we display the likelihood contours in the parameter space in the top left and bottom two plots of Figure 6 for a bin centered at r = 16.93 kpc, where the maximum dip in the β profile was seen in Figure 5. Additionally, in the corresponding bin, the top right plot in the figure demonstrates the posterior distributions of the model parameters σ_r, σ_θ, and σ_ϕ of the GVE model obtained from 5 × 10⁵ MCMC random walks. One can see that even at a distance of just twice that of R_☉, the σ_θ and σ_ϕ distributions are quite broad when compared to the σ_r distribution; this is the reason for the large uncertainty in the value of β as we move outward in r.

Zoom In Zoom Out Reset image size

3.4. Effect of the Halo Substructures

There is now enough observational evidence to suggest that the stellar halo is highly structured, particularly as one moves outward into the halo (Bell et al. 2008). Using clustering algorithms on simulated N-body stellar halos, Sharma et al. (2011b) find that the fraction of material in substructures increases monotonically as a function of distance from the center and at around 65 kpc can be as high as 50%. Hence, should we include the substructures or exclude them while studying the kinematic properties of the halo? If the kinematic properties of a sample are dominated by a few massive accretion events, then one should exclude them. However, in spite of being highly structured, if the sample is a superposition of a large number of events with none of them being individually too dominant, then it is best to include them. Results of Sharma et al. (2011b) on simulated halos show that for the range of radii that spans our BHB stars (r < 40 kpc), the amount of material in the substructure should be less than 20%. So we expect the substructures to have a minor effect on the kinematic properties that we have derived. However, it is still important to check if this is true.

To study the effect of substructures on our estimation of the dispersions, we mask two prominent features of the halo that contaminate our sample, namely, the Sagittarius stellar stream and the Virgo overdensity. The cuts we impose include the Lambert equal-area projection cut as given in Bell et al. (2008) and an additional cut in equatorial coordinates (RA and Dec.) suggested by Deason et al. (2011b). We mask the region within 0 < X (abscissa of the equal-area projection) <30, where X is given by 63.63961[2(1–sin b)]^1/2; and 0° < RA < 50°, and −30° < Dec. < 0°, which is a purely geometric cut. These stringent cuts reduce our final sample to 2975 stars. A proper phase-space masking of these structures will be revised in future work.

In Figure 7, we present our results obtained after masking the substructures. As masking reduces the sample size almost by half, we employ CME with n_bin = 500 stars only, instead of n_bin = 1200 as was done earlier with the unmasked data, to avoid excessive smoothing of the estimated profiles. Figure 7 shows that the masking of substructures has little effect on the estimation of velocity dispersion profiles, and the β profile is almost unchanged. This alleviates the concern that perhaps the turnover points in the velocity dispersion profiles and the dip in the β profile, seen in Figure 3, in the region r = 13–18 kpc could be due to the dominance of halo substructures.

Zoom In Zoom Out Reset image size

In their recent works, D11 and D12 fitted a DF of the form given by Equation (3) to the BHB samples obtained from the SEGUE survey in order to estimate the model parameters. The models adopted are constant anisotropy models given by Equation (2). In D11, the potential is assumed to be a power law (∝r^−γ) with a constant index (γ = 0.35). Later in D12, they break the degeneracy in their model and consider the potential normalization (Φ₀), potential slope (γ), and anisotropy (β) as free parameters. Note that in D11, there is an additional parameter specifying rotation (odd part of the DF), but this was dropped in D12 as they were not focusing on rotation. The methodology applied to measure the model parameters is similar to ours, which involves marginalizing the DF over the tangential velocities to derive the LOSVD: fitting the LOSVD to the data using the maximum likelihood method and in return obtaining the best estimates of the model parameters.

4.1. Anisotropies of Inner and Outer Halos by D11

4.2. Anisotropy at 16 < r/kpc < 48 Seen by D12

D12 re-calibrate the distances of BHB stars using the color–magnitude relation given by Deason et al. (2011b) and then select stars within 16 ≲ r/kpc ≲ 48 (1933 stars) from the X11 BHB samples. They fit the D11 model to study the nature of β in the outer halo. Unlike D11, as mentioned previously, here they break the degeneracy found in their model and consider β, Φ₀, and γ as free parameters. Thus, while they fit the model, they simultaneously estimate all three parameters. They find β = 0.4^+0.2_−0.2 for α = 4.6. Using a model allowing for oblateness (q = 0.59), they find β = 0.5^+0.1_−0.2. If we apply the GVE DF to stars in the range 16 < r/kpc < 48, we find β = −0.14^+0.52_−0.66 (σ_r = 97.3^+2.9_−3.0 km s⁻¹, σ_θ = 122.7^+26.4_−33.0 km s⁻¹, and σ_ϕ = 78.5^+34.3_−40.7 km s⁻¹, posterior distributions are shown in Figure 8). This more or less looks like the mean value of β in this range (see Figure 5), provided we take into account the fact that the number density of stars peaks at around r = 16 kpc. Although the D12 value is still within our 1σ region, our predicted value is lower than that in D12. Figure 5 shows that β is not constant in the range 16 < r/kpc < 23. It increases from being tangential to isotropic. Using the D11 DF also gives similar results (Section 3.3). Beyond this range, we cannot directly measure β, but by deriving a best-fit circular velocity profile and making use of the σ_r profile which is available until r = 56 kpc, we can predict β, and we return to this in Section 5. However, an assumption about α also has to be made. Beyond r > 27 kpc, the density slope has been shown to change from −2.4 to −4.5. Adopting a steeper density slope increases the value of β. For α = 4.5, we find β ∼ 0.5 for r > 27 kpc; this is more in agreement with D12. To conclude, the D12 value of β = 0.4 for 16 < r/kpc < 48, although derived for a sample which is dominated by stars within r < 27 kpc, is not appropriate for the range 18 < r/kpc < 23, instead it is correct for the range 23 < r/kpc < 48.

Zoom In Zoom Out Reset image size

Finally, we check how we can best measure β in the outermost parts, 35 < r/kpc < 84, using the D11 DF and the data in hand. This region consists of 762 stars, and we assume α = 4.5 and γ = 0.35 and repeat the analysis as in D12. The likelihood distributions of model parameters are shown in Figure 9. The mass–anisotropy degeneracy is clearly visible here, suggesting that it is very difficult to directly measure β unless an explicit or implicit assumption about the potential is made.

Zoom In Zoom Out Reset image size

Here, we use the measured values of σ_r(r) and β(r) from our analysis given in Figure 3 to determine the circular velocity curve of the Galaxy (v_circ) through the Jeans equation (Equation (7)). Besides anisotropy and radial velocity dispersion information, we also need to adopt some density profile for the tracer population but not of the spectroscopic sample. We adopt a value of α = 2.4, as suggested by recent works of Watkins et al. (2009) and Deason et al. (2011b), for the range of distance explored here (r < 25 kpc). The blue dots with error bars in Figure 10 are our estimates of v_circ using CME with n_bin = 1200. The uncertainties on v_circ were computed using a Monte Carlo-based scheme from uncertainties in σ_r and β. The leftmost and rightmost points have comparatively larger error bars than intermediate points. For the leftmost point, the large error bar is mainly due to the large error in the value of σ_r. On the other hand, for the rightmost point, the large error bar is mainly due to large error in the values of σ_θ and σ_ϕ.

Zoom In Zoom Out Reset image size

Figure 10 shows that the circular velocity profile derived from our analysis (blue points) displays prominent features. We now check if such features are also observed in other studies using tracers other than BHB stars. For this, we overplot the v_circ compiled by Sofue et al. (2009) as black dots, obtained from several references (Burton & Gordon 1978; Blitz et al. 1982; Clemens 1985; Fich et al. 1989; Honma & Sofue 1997a, 1997b; Honma et al. 2007; Demers & Battinelli 2007). Further details about the source of each individual point can be found in Sofue et al. (2009) and references therein. Note that the v_circ values of Sofue et al. (2009) are computed for (R_☉, v_LSR) = (8.0 kpc, 200.0 km s⁻¹). Correcting v_circ for our adopted values of (R_☉, v_LSR) = (8.5 kpc, 220.0 km s⁻¹) is beyond the scope of this work and, hence, we simply overplot these published values in Figure 10. One can see that there is a prominent dip at 9 kpc in the Sofue et al. (2009) compiled v_circ profile. They explain this dip by introducing massless rings on top of a disk with exponential surface density. We also find a similar dip in our v_circ profile, but at around 10–12 kpc. The slight shift in the position of the dip could be due to the large width of our bins and also the fact that unlike Sofue et al. (2009), who measure v_circ in the midplane of the Galaxy, we measure v_circ over a spherical shell that intersects with the SDSS footprint.

We now estimate the mass of the dark matter halo of the Milky Way, assuming a three-component model of the Galaxy consisting of the bulge, the disk, and the halo. The bulge is modeled as a Hernquist sphere and the disk is assumed to follow an exponential profile (Xue et al. 2008). The parameters for the bulge and disk are taken from Sofue et al. (2009) and are kept fixed. Although Sofue et al. (2009) use massless rings, here we have ignored them since our main aim is to fit the dark matter halo. We model the dark matter halo using the Navarro–Frenk–White (NFW) density profile (Navarro et al. 1996). Here, we consider both the halo and the bulge to be spherically symmetric. The non-axisymmetric effect due to a bar-shaped bulge is neglected here. Potentials for the bulge, the halo, and the exponential disk can be expressed, respectively, as

$$\begin{eqnarray} \Phi _{\rm disk} (r)& = &{-}\frac{GM_{\rm disk}(1-e^{-r/b})}{r} , \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \Phi _{\rm bulge} (r)& = &{-}\frac{GM_{\rm bulge}}{r+a} ,{\rm } \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} \Phi _{\rm NFW} (r) & = &{-}\frac{GM_{\rm vir}\ln (1+r c/R_{\rm vir})}{g(c)r}, \end{eqnarray} \tag{ 12 }$$

where M_disk = 6.5 × 10¹⁰ M_☉, b = 3.5 kpc, M_bulge = 1.80 × 10¹⁰ M_☉, and a = 0.5 kpc (Sofue et al. 2009).

Note that the disk potential as given by Equation (10) is spherically symmetric. It means that the disk is considered to be a spherical body with an exponential surface density falloff. To get an idea of the error that is incurred due to the assumption of the disk as a spherical body with the same mass as the flattened disk, we refer the reader to Binney & Tremaine (2008, Figure 2.17). The maximum error in v_circ is roughly 13%, which is at a distance of about twice that of the disk scale length. However, at the larger distances along the midplane, the discrepancy is smaller. In the general case of triplanar symmetry (elliptic disk), in reality, the disk potential has to be the function of both polar and azimuthal coordinates, and in the special case of axial symmetry (circular disk) it has to be the function of sole polar coordinates, in addition to the radial coordinates (r). We use the spherically symmetric form here for two reasons. First, we make use of the spherical form of the Jeans equation given by Equation (7), which demands a spherical potential. Second, it is to ease the comparison with earlier studies, e.g., Xue et al. (2008), that adopt a similar definition. However, later on, we consider a three-dimensional (3D) disk potential and discuss its consequences for the estimation of mass. The function g in the NFW potential is given by

$$g(c) = \ln (1+c) - \frac{c}{1+c}$$

and

$$R_{\rm vir} =\left(\frac{2 M_{\rm vir}G}{H_0^2 \Omega _m \Delta _{\rm th}}\right)^{1/3}.$$

The total potential Φ(r) of the Galaxy is then simply

$$\begin{equation} \Phi (r) = \Phi _{\rm bulge}(r) + \Phi _{\rm disk}(r) + \Phi _{\rm NFW}(r). \end{equation} \tag{ 13 }$$

We adopt the value of the Hubble constant, H₀ = 70.4 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 (Komatsu et al. 2011), and Δ_th = 340 (Bryan & Norman 1998).

A NFW halo has two free parameters, the mass M_vir and the concentration c. Since we do not have enough data points spanning a wide range in radius, we avoid fitting both the parameters simultaneously. Instead, we use the concentration–mass relation

$$\begin{equation} c = 327.3\ M_{\rm vir}^{-0.12}, \ \ 10^{11} \leqslant M_{\rm vir}/M_{\odot } \leqslant 10^{13}, \end{equation} \tag{ 14 }$$

as has been reported in N-body simulations of dark matter halos (Macciò et al. 2007).

Finally, we can derive the resultant circular velocity (v_circ) from the total potential (Equation (13)) by computing (rdΦ/dr)^1/2. We fit the obtained theoretical rotation velocity curve to our observed values of v_circ shown by the blue dots in Figure 10, and the red line is our best-fit circular velocity curve. The v_circ profiles for the different components are also shown separately. The dashed black line is the corresponding best-fit NFW halo profile. The best-fit value for the fit parameter, M_vir, for our three-component baryon and dark matter mass distribution is 0.9^+0.4_−0.3 × 10¹² M_☉. The corresponding values of R_vir and c derived from the best-fit value of M_vir are 249⁺³⁴₋₃₁ kpc and 12.0^+0.6_−0.5, respectively. We estimate the mass of the Galaxy within r ≲ 25 kpc to be 2.1 × 10¹¹ M_☉. Assuming a functional form for P(v_los/v_circ) obtained from simulations, Xue et al. (2008) derive v_circ from the v_los of BHB stars. The derived v_circ is then used to estimate the virial mass of the dark matter halo. They find M_vir = 0.91^+0.27_−0.18 × 10¹² M_☉, which is in good agreement with our result; note that uncertainties are, however, slightly larger in our estimates. Unlike them, we do not make any assumption about the functional form of P(v_los/v_circ).

Here, we study the effect of chosen disk models for which we consider a more realistic 3D potential for the disk by Miyamoto & Nagai (1975), which holds for the special case of a circular disk. The expression for this potential is given by

$$\begin{equation} \Phi _{\rm disk} (R,z) = -\frac{G M_{\rm disk}}{\sqrt{R^2+(a+\sqrt{z^2+b^2})^2}}. \end{equation} \tag{ 15 }$$

Here, again, the disk parameters are obtained from the best-fit values along the galactic midplane (z = 0), which reproduces the v_circ(R) profile for the Sofue et al. (2009) razor-thin exponential disk (Freeman 1970), i.e., b = 0.0 kpc.² Hence, the best-fit value for the disk parameter is a = 2.5 kpc, whereas the mass M_disk = 6.5 × 10¹⁰ M_☉ is the same as in Sofue et al. (2009). Since we assume the disk potential to be three-dimensional here, for the purpose of computing the total v_circ, we consider the component of the force along the radial direction (r) only. The bulge and the halo models are the same as earlier. The best-fit values for the NFW halo parameters with the updated disk model are found to be M_vir = 1.2^+0.5_−0.4 × 10¹² M_☉, with R_vir = 274⁺³⁵₋₃₀ kpc. Instead, if we consider the total magnitude of the force in order to compute v_circ, then we estimate M_vir(R_vir = 269⁺³⁴₋₃₂ kpc) = 1.1^+0.5_−0.4 × 10¹² M_☉. Historically, large values of a have been assumed while modeling the disk (Wolfire et al. 1995; Bland-Hawthorn 2009). We find that this leads to a much more massive dark matter halo—for a = 4.5 kpc we find M_vir(R_vir = 299⁺³⁶₋₃₃ kpc) = 1.6^+0.6_−0.5 × 10¹² M_☉ and for a = 6.5 kpc we find M_vir(R_vir = 321⁺³⁵₋₃₄ kpc) = 1.9^+0.7_−0.6 × 10¹² M_☉.

There are few things that have insignificant or unexplored effects on our mass estimation, e.g., we do not take into account the mass of the supermassive black hole, which is ∼4 × 10⁶ M_☉ (Schödel et al. 2002; Eisenhauer et al. 2005; Gillessen et al. 2009), and is approximately 1/1000 of the mass of the bulge. Moreover, its effect is like that of a point mass and can be easily absorbed into the bulge mass. Another effect that is not considered is the tidal effect on the primary object that has been qualitatively studied to find the mass ratio between the Galaxy and M31 (Baiesi Pillastrini 2009) and depends strongly on different impact parameters (Eneev et al. 1973). We find that the R_vir of the Galaxy is ∼250 kpc and for M31 it is ∼260 kpc (Seigar et al. 2008; Majewski et al. 2007). Given that the distance between these two galaxies is ∼780 kpc (Ribas et al. 2005; McConnachie et al. 2005; Karachentsev & Kashibadze 2006), which is more than double the sum of their virial radii, we believe that the tidal effect of M31 on the overall mass estimate of the Galaxy, if any, should be negligible. The tidal effects of the LMC and the SMC on the Galaxy have not been explored in this paper.

Note that our mass modeling of the Milky Way does not make any assumption about the value of β, instead we use the value of β directly computed from the data. The only assumption that we make is that the density of the dark matter halo follows an NFW profile. As long as that assumption holds, our estimates for v_circ and the mass of the Milky Way should also be valid in the outer parts r > 25 kpc, where we cannot directly measure β. If one wants to directly measure v_circ in the outer parts using only line-of-sight velocities, then one has to make an assumption about the underlying β. Several attempts have been made in this regard, with each of them making different assumptions about β, and hence introducing a bias into the estimated mass. Below, we compare these with our prediction for the mass of the Galaxy. The dashed line in Figure 11 is our best-fit model of v_circ. Since literature sources mostly report mass within a certain radius, to facilitate comparison we convert it to v_circ using the relation

$$M(<r) = v_{\rm circ}^{2}r/{\rm G}.$$

In Figure 11, the plotted v_circ from different sources span a wide range in radii and were computed using different types of tracer populations. By fitting a model to the kinematics of the satellite galaxies and globular clusters, Wilkinson & Evans (1999) measure the mass to be M(50 kpc) ∼ 5.4^+0.2_−3.6 × 10¹¹ M_☉. This agrees with estimates by Kochanek (1996) (M(50 kpc) = (4.9 ± 1.1) × 10¹¹ M_☉) and Sakamoto et al. (2003) (M(50 kpc) = 1.8–2.5 × 10¹¹ M_☉). Watkins et al. (2010) apply the tracer mass estimator formalism to 26 satellite galaxies and find that M(300 kpc) = (0.9 ± 0.3) × 10¹² M_☉. Their mass estimate is, however, prone to the systematics introduced by using an assumed β, as duly mentioned by them. One can see in the figure that at r = 100 kpc, depending upon the chosen anisotropy, their mass could vary anywhere between 0.3 × 10¹² and 1.4 × 10¹² M_☉. Studying the hypervelocity stars within 80 kpc, and assuming β = 0.4, Gnedin et al. (2010) estimate M = 6.9^−1.2_+3.0 × 10¹¹ M_☉, which is slightly higher than our estimate. Using BHB stars and the Watkins et al. (2010) tracer formalism, Samurović & Lalović (2011) estimate the mass at r = 85 kpc to be 8.83 ± 0.73 × 10¹¹ M_☉. This is slightly higher than our estimate, probably because they assume β = 0. With the mixed sample of tracers (BHB and CN stars) populating the outermost halo (r ∼ 50–150 kpc), Deason et al. (2012b) estimated the mass of the Galaxy to be M(150 kpc) = (5–10) × 10¹¹ M_☉. The variation is mainly due to the uncertainty related to the adopted potential and density slopes and anisotropy. Their range of mass at the outermost halo falls within our estimation.

Zoom In Zoom Out Reset image size

We now study the kinematic properties of simulated stellar halos in which halos are formed purely by accretion of satellite galaxies. For this, we use the simulations of Bullock & Johnston (2005). In order to construct a synthetic sample of BHB stars from these simulations, we use the code GALAXIA (Sharma et al. 2011a). Figure 12 shows the velocity dispersion and anisotropy profiles of 11 different ΛCDM halos as a function of galactocentric radius r. The mean of all the halos is also shown alongside the thick red line. In general, the velocity dispersions fall off with radius. At small r, the fall is rapid, but at large r it is much slower. Asymptotically, the ratio σ_r/V_vir approaches a value of around 0.8. Note that β rises rapidly from a value of zero in the center to about 0.8 at r ∼ 10 kpc and, thereafter, shows very little change. These results are in good agreement with results of Abadi et al. (2006) and Sales et al. (2007), who study stellar halos formed in cosmological hydrodynamical simulations including star formation and feedback.

Zoom In Zoom Out Reset image size

First, note that ΛCDM halos are rarely tangential for any given radius. For most of the range of r, β is in general greater than 0.5. We fit an analytic function of the Cuddeford (1991) form given by

$$\begin{equation} \beta (r)=\beta _0 \frac{r^2}{r^2+r_0^2} \end{equation} \tag{ 16 }$$

to the mean β profile of the 11 ΛCDM halos. The best-fit values of the free parameters were found to be β₀ = 0.765 and r₀ = 0.00843 R_vir. Figure 13 shows that the fit is quite good for a wide range of r. The slight mismatch at r < 1 kpc could be due to issues related to force resolution. In the outer parts, most of the mass is in bound structures and hence is not smoothly distributed. This is probably responsible for the non-smooth behavior in β in the outer parts.

Zoom In Zoom Out Reset image size

Figure 14 presents the beta profiles for halos with a non-ΛCDM accretion history, i.e., halos having accretion history significantly different from that predicted by the ΛCDM model of galaxy formation. Six halos that we consider have accretion events that are dominated by (1) radial orbits, (2) circular orbits, (3) old events, (4) recent events, (5) higher luminosity, and (6) low luminosity, and these are from simulations of Johnston et al. (2008). Signatures of different accretion events can be seen in the β profiles. The most significant difference is between the radial and the circular halo, which is expected since the orbital properties of the satellites were different to begin with. Note that the circular halo is the only one among all the simulated halos that can have β < 0. The old halo almost perfectly follows the mean profile that we had derived for the 11 ΛCDM halos and has the smoothest profile. This is due to the fact that the stars in this halo are completely phase mixed and have no structures of any kind. The young halo has very few stars in the inner regions and shows non-smooth behavior due to the presence of a significant number of structures. The high-luminosity halo is also very similar to the old halo. However, the low-luminosity halo has higher β for r < 0.01 R_vir. This is most likely due to circularization of orbits when acted upon by dynamical friction. Orbit circularization has also been reported by Sales et al. (2007) in their simulations. Note that the effect of dynamical friction is strongest for high-luminosity events and weakest for low-luminosity events. Moreover, when acted upon by dynamical friction, satellites lose energy and move toward the inner regions of the halo. This partly explains why the high-luminosity halo has low β in the inner regions as compared to the low-luminosity halo.

Zoom In Zoom Out Reset image size

In Section 3.2, we measured β until r = 23 kpc. Beyond this, we can measure σ_r out to r = 56 kpc, but not σ_ϕ and σ_θ. By using the circular velocity curve that we derived in Section 5, we can predict β beyond r > 23 kpc by making use of the Jeans equation (Equation (7)). To proceed, we need to make an assumption about the density slope (α); beyond r > 27 kpc, it has been shown that α is around 4.5. Assuming this value, the predicted β is plotted in Figure 15. One can see that there is a slight jump in the value of β passing from r = 23 to 27 kpc, and beyond this the value of β is around 0.4. The sudden jump in β occurs via the Jeans equation (Equation (7)), due to the discontinuity in α (=−dln ρ/dln r). Note that an assumption of the steeper density slope would increase β and vice versa. The red line in the figure shows the anisotropy profile fitted to the simulated ΛCDM halos given by Equation (16) from Section 6. It can be seen that accretion-based models cannot explain the dip in the observations, especially the profile in the region 12 < r/kpc < 23. However, outside this region, the simulations are roughly in agreement with observations. For r < 12 kpc, the observations match the value of around β = 0.5 seen in simulations. For r > 23 kpc, the overall value of β in observations is slightly low but the profile is flat as in the simulations. The outer halo at r = 56 kpc is radial with β = 0.55.

Zoom In Zoom Out Reset image size

We study the kinematics of ∼4500 BHB stars to obtain the velocity dispersion profiles along three orthogonal axes in spherical polar coordinates using the GVE DF. GVE as an estimator of the velocity dispersion has the advantage that no assumptions about potential or density are needed a priori. From the estimated velocity dispersion profiles using maximum likelihood analysis, we also derive the anisotropy profile of the Galactic halo and compare it to the simulated ΛCDM halos. Finally, using the radial velocity dispersion profile, anisotropy profile, and density power law, we constrain the mass of the Milky Way using the Jeans equation.

We measure the σ_r profile of the halo out to r ∼ 60 kpc. At large distance (d ≫ R_☉), σ_r can be approximated by σ_los. Thus, in the outskirts, the σ_los profile given by Xue et al. (2008) converges to our σ_r profile. At r ∼ 60 kpc, σ_r attains ⋍100 km s⁻¹. However, in the inner halo, the approximation (σ_r ≈ σ_los) is invalid and we find that the deviation of σ_los from σ_r is as high as ∼40 km s⁻¹. We obtain a σ_r(r) profile with a plateau in the inner halo, a sudden fall at r ∼ 15 kpc, and a gradual decline outward. Qualitatively, a similar profile is also found by Sommer-Larsen et al. (1997). However, our σ_r profile sharply falls at r ∼ 15 kpc, whereas they find a gentle transition. This is probably due to the fact that they make an assumption that v_circ(r) is constant, which we find is not completely true.

Next, we estimate the tangential velocity dispersions, σ_θ and σ_ϕ. Using these estimates, we are able to measure β(r) until r = 25 kpc. Astonishingly, we discover a dip in the β profile at r = 17 kpc, where β⋍ − 1.2. We find that the inner halo (r < 12 kpc) is radial with β⋍0.5. This result of the radially biased inner halo concurs with the recent results by Smith et al. (2009a) and Bond et al. (2010) using the proper motions. Beyond the switchover point in the range 18 ≲ r/kpc ≲ 25, the anisotropy rises slightly and becomes isotropic to mild radial. We also verify the result using an alternative DF, namely, the D11 DF. A small systematic in the β(r) profiles is seen from these two models, which is mainly due to the assumption about the density slope (α) that needs to be made a priori for D11 DF. We check for the contribution of the halo substructures, namely, the Virgo overdensity and Sagittarius stellar stream, and find that they have little effect on the anisotropy profile. The effects of v_LSR and R_☉ upon our velocity dispersions and anisotropy estimates are also found to be negligible.

D11 study the BHB stars in the radial bin (10 < r/kpc < 25) and find the halo to be tangential. After reanalyzing the D11 sample in this bin, we find that this is mainly because of the choice of their bin size that encompasses the transition region (13 < r/kpc < 17) where we detected a dip in the β(r) profile. Possibly, it could also be because of the degeneracy between the potential and anisotropy in their model. However, in their recent work, D12 break the degeneracy among their model parameters and measure β = 0.4^+0.15_−0.2 with α = 4.6 in the region 16 ≲ r/kpc ≲ 48. Within the range of uncertainty, our value for β (−0.14^+0.52_−0.66) using the GVE model agrees with them. We find that the D12 value of β ≈ 0.4 in this bin, although derived from a sample which is dominated by stars within r < 27 kpc, is not appropriate for the range 18 < r/kpc < 23, instead it is appropriate for the range 23 < r/kpc < 48. Finally, we check how well we can estimate the mass and anisotropy together using D11 DF in the outermost region (35 < r/kpc < 84). We find that due to the lack of tangential information, a degeneracy between mass and anisotropy cannot be broken.

Substituting the estimates of σ_r(r) and β(r) in the Jeans equation, we then calculate the circular velocity profile of the Galaxy (v_circ(r)). We detect the dip in the v_circ profile at 10–12 kpc, also seen by Sofue et al. (2009) at 9 kpc, which is attributed to the massless ring as a perturbation to the disk. Finally, we fit the three-component (exponential disk, Hernquist bulge, and NFW halo) galaxy model to the observed v_circ profile in order to obtain the mass distribution of the Galaxy. From our best-fit model, we calculate M_vir of the halo to be 0.9^+0.4_−0.3 × 10¹² M_☉ with R_vir = 249⁺³⁴₋₃₁ kpc and the concentration parameter, c = 12.0^+0.6_−0.5. The mass of the Galaxy, within the extent that we are able to constrain all the three components of velocity dispersions, is estimated to be M(r ≲ 25 kpc) = 2.1 × 10¹¹ M_☉. Our estimate for M_vir is in good agreement with most of the recent estimates, namely by Watkins et al. (2010), Kochanek (1996), Wilkinson & Evans (1999), Sakamoto et al. (2003), Sofue et al. (2009), and Deason et al. (2012b), as demonstrated in Figure 11. In their studies of the same population of stars (BHB), Xue et al. (2008) also fit a three-component galaxy model and calculate M_vir to be 0.91^+0.27_−0.18 × 10¹² M_☉. Our result for M_vir is in very good agreement with their estimate but our uncertainty is slightly larger. Note that they make an assumption about (v_circ/v_los) with radius from the simulations, whereas we do not make any such assumption. Additionally, we also consider a more realistic 3D disk model that is found to predict a slightly higher Galactic mass M_vir = 1.2^+0.5_−0.4 × 10¹² M_☉ with R_vir = 274⁺³⁵₋₃₀ kpc for flattening constant a = 2.5.

In the end, we used the measured quantities σ_r and v_circ to extend the β profile beyond r ∼ 25 kpc up to the distance where σ_r(r) can be confidently measured (r ∼ 60 kpc). The only assumption that we make here is about the density profile which we set at α = 4.5 in consonance with the recent results by Watkins et al. (2009) and Deason et al. (2011b). We find that the outer halo is radial and attains β = 0.55 at r ∼ 60 kpc.

We also compare our result with the simulated stellar halo, which is formed purely by accretion (Bullock & Johnston 2005). This simulated halo is found to be in rough agreement with the observed halo in the inner region r < 12 kpc. One can see that none of the simulations the β profiles obtained could predict a tangential halo at any distance, and thus fails to explain a dip seen at r = 17 kpc in the observed β profile. In contrast, in the outer region (r > 25 kpc), simulations and observations both agree in an overall sense with the anisotropy and predict a flat anisotropy profile.

In all of the observed quantities σ_r, σ_θ, σ_ϕ, and β, we see a dramatic shift in properties at r ∼ 17 kpc. We noted that these undulations in the profiles are translated into our v_circ estimation, resulting in a varying v_circ profile. It could be true the other way around, in the sense that the non-monotonic trends seen in all of our kinematic profiles could be due to the presence of so far unaccounted features in the Milky Way potential.

Alternatively, the shift in the properties seen in the observed profiles could possibly be an indication of a complex multi-component halo. Recently, there have been series of works advocating a multi-component halo. The studies of the calibration stars by Carollo et al. (2007, 2010) from the SDSS survey and the follow-up studies by Beers et al. (2012) have shown that the halo has at least two distinct components. They determined the inner halos to be formed in situ, whereas the outer halos are considered to be formed by accretion. Carollo et al. (2010) and de Jong et al. (2010) have found that the population fraction inversion point between the inner and outer halos lies between ∼15 and 20 kpc. Kinman et al. (2012) studied the BHB and RR Lyrae stars toward the Galactic anti-center and North Galactic Pole and found that the retrograde component of the halo dominates for r > 12.5 kpc. It seems that this transition between the inner to the outer halos is recorded in the β of the halo as well. Additionally, duality in the formation history of the halo has also been seen in the recent smooth particle hydrodynamics and N-body simulations by Zolotov et al. (2009), McCarthy et al. (2012), and Font et al. (2011). Kalirai (2012) recently attributed the age difference of two billion years in the halo components to in situ accretion. On the contrary, Schönrich et al. (2011) reanalyzed the calibration stars from Carollo et al. (2010) and found no reliable evidence for the existence of an outer retrograde halo.

In a nutshell, the stellar halo is a test bed to understand the formation history of the galaxy. Even in this era, where we have access to a huge volume of spectroscopic and photometric data, the crucial physical quantities like velocity dispersions and anisotropy are not completely understood due to the lack of proper motions. With the advent of data in the coming decades from the magnificent next generation of spectroscopic surveys like LEGUE (Deng et al. 2012) and, especially, unprecedented proper motions from an astrometric mission like Gaia (Perryman 2002), we will be able to put strong constrains on these fundamental quantities. Additionally, in order to see the bigger picture, we will need to confirm the results with different stellar types or an alternative tracer. Moreover, exploring the southern sky is equally important in order to complete the picture, for which upcoming spectroscopic surveys like GALAH³ will also play an important role.

We sincerely thank the anonymous referee for comments that helped to improve the paper. We also thank Dr. Ralph Schonrich and Francesco Fermani for their comments on the manuscript, particularly on the effect of (v_LSR, R_☉). Dr. Xiang Xiang Xue is thanked for the online publication of the clean sample of BHB stars. Tim White is also thanked for comments on the original manuscript. P.R.K. acknowledges the University of Sydney International Scholarship (USydIS) for the support of his candidature. G.F.L. acknowledges support from ARC Discovery Project DP0665574. J.B.-H. is funded through a Federation Fellowship from the Australian Research Council (ARC) and S.S. is funded through ARC DP grant 0988751 which supports the HERMES project.

Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web site is http://www.sdss.org/.

Here, we investigate an effect of n_bin on our analysis. Figure 16 shows the velocity dispersion profiles and the anisotropy profile for the same sample of stars but with different particles in each bin. We see that with the decrease in the number of stars in each bin, the uncertainties in the result increase. However, the overall trend remains unaltered.

Zoom In Zoom Out Reset image size

In the literature, there are varied claims about the value of v_LSR ranging from 184 to 270 km s⁻¹ (Olling & Merrifield 1998; Méndez et al. 1999; Reid et al. 2009; Bovy et al. 2009; Koposov et al. 2010). Similarly, the value of R_☉ is also debatable within 8–8.5 kpc (Reid 1993; Ghez et al. 2008; Gillessen et al. 2009). McMillan & Binney (2010) found that the ratio of v_LSR/R_☉ can be better constrained than each of them alone and should range between 29.9 and 31.6 km s⁻¹ kpc⁻¹. Distressingly, there is still no consensus upon the values of (v_LSR, R_☉). To study the effect of chosen values of (v_LSR, R_☉) upon our estimates of dispersion profiles, we repeat the same analysis performed to obtain the black diamond points in Section 3 (Figure 3). In Figure 17, we show the results for different values of (v_LSR, R_☉). Here, again, the black diamond markers are obtained for a case (v_LSR, R_☉) = (220.0 km s⁻¹, 8.5 kpc), and the figure is thus a replica of the diamond points from Figure 3, repeated for ease of comparison. The red and black markers in the figure show the effect of chosen R_☉ values upon our estimates, whereas the red, blue, and cyan markers demonstrate the effect of chosen v_LSR. Within the range of (v_LSR, R_☉) investigated, the figure clearly demonstrates a negligible effect upon our estimates, given the range of uncertainties.

Zoom In Zoom Out Reset image size