1 Introduction

On the Shoulders of Giants: Properties of the Stellar Halo and the Milky Way Mass Distribution


Halo stars orbit within the potential of the Milky Way and hence their kinematics can be used to understand the underlying mass distribution. However, the inferred mass distribution depends sensitively upon assumptions made on the density and the velocity anisotropy profiles of the tracer population. Also, there is a degeneracy between the parameters of the halo and that of the disk or bulge. Most previous attempts that use halo stars have made arbitrary assumptions about these. In this paper, we decompose the Galaxy into 3 major components – a bulge, a Miyamoto-Nagai disk and an NFW dark matter halo and then model the kinematic data of the halo Blue Horizontal Branch and K-giant stars from the Sloan Extension for Galactic Understanding and Exploration (SEGUE). Additionally, we use the gas terminal velocity curve and the Sgr A proper motion. With the distance of the Sun from the centre of Galaxy , our kinematic analysis reveals that the density of the stellar halo has a break at , and an exponential cut-off in the outer parts starting at . Also, we find the tracer velocity anisotropy is radially biased with in the outer halo. We measure halo virial mass to be , concentration to be , disk mass to be , disk scale length to be and bulge mass to be . The mass of halo is found to be small and this has important consequences. The giant stars reveal that the outermost halo stars have low velocity dispersion but interestingly this suggests a truncation of the stellar halo density rather than a small overall mass of the Galaxy. Our estimates of local escape velocity and dark matter density ( GeV cm) are in good agreement with recent estimates. Some of the above estimates, in particular , are depended on the adopted value of and also, on the choice of the outer power-law index of the tracer number density.

Subject headings:
galaxies: individual (Milky Way)-Galaxy: halo - stars: giants - stars: kinematics

1. Introduction

Mass is the fundamental property of any galaxy. An accurate measurement of the Galaxy mass has repercussions in many sectors, e.g., in its mass assembly history (Wechsler et al., 2002), identifying a realistic Galaxy in a simulation(e.g., Vera-Ciro et al., 2013; Wang et al., 2012) or its analogue (Robotham et al., 2012), simulating the tidal streams or the orbit of the satellite galaxies (e.g., Newberg et al., 2010), studying the tidal impact of the Galaxy on the satellite galaxies (e.g., Johnston, 1998; Kallivayalil et al., 2013; Nichols et al., 2014) etc. Various approaches have been undertaken to determine the mass distribution of the Galaxy, e.g., the timing argument (Kahn & Woltjer, 1959; Li & White, 2008), the local escape speed (Leonard & Tremaine, 1990; Smith et al., 2007; Piffl et al., 2014), the orbital evolution of the satellite galaxies and globular clusters (Lin & Lynden-Bell, 1982; Boylan-Kolchin et al., 2013), modeling the tidal streams (Law et al., 2005; Newberg et al., 2010; Sanders & Binney, 2013b; Sanderson & Helmi, 2013), the HI gas rotation curve (e.g. Sofue et al., 2009), fitting a parametrized model to the available observational constraints (e.g. Dehnen & Binney, 1998; McMillan, 2011; Irrgang et al., 2013) etc. Each of this method has its own inherent limitation, for example, the local escape speed method suffers from the paucity of high velocity stars and also it is unclear whether the phase space is fully filled up to the escape velocity. The use of HI gas rotation curve suffers from the fact that there is no extended HI disk reported for the Galaxy and hence fails to probe the mass that lies beyond the extent of the disk. For an in-depth discussion and description of the various methods we refer the reader to the two recent reviews: by Courteau et al. (2013) on the galaxy masses and by Sofue (2013) on the rotation curve and references therein.

A simple, yet robust method to probe the Galaxy mass is provided by an application of the Jeans (1915) Equation. The spherical Jeans Equation for a system in dynamic equilibrium is given by


where is gravitational potential, is stellar density, is radial velocity dispersion and is velocity anisotropy, and all of them can be a function of galactocentric distance . Thanks to massive stellar surveys such as SDSS/SEGUE that provide a catalog of position and radial velocity measurements of large number of halo tracers, it is now possible to use the Jeans Analysis to put the most stringent constraints on the halo parameters out to the maximum observed distance. Among all the parameters that enter the Jeans Analysis, the most uncertain quantity is velocity anisotropy , which is defined as


where and are the velocity dispersions along the spherical polar () and azimuthal () directions. With , signifies dominance of radial motion of stars, signifies dominance of tangential motion and an isotropic system. In is not known then the Jeans Analysis suffers from a degeneracy known as the mass-anisotropy degeneracy. In general terms it means that same radial velocity dispersion profile can be obtained either by lowering the value or by increasing the mass. Numerous studies based on the kinematics of different stellar species, namely the sub dwarfs (Smith et al., 2009b), the main-sequence stars (Brown et al., 2010) and the BHB stars (Kafle et al., 2012, hereafter K12) concur that (radial) in the Solar neighborhood. There have been recent attempts to constrain the beyond the Solar neighborhood, including K12, who use the line-of-sight velocity of a BHB sample to measure to a radius of . They find a non-monotonic trend in starting with 0.5 (radial) at small which falls to -1.2 (tangential) at and then rises again to 0 at . An additional measurement of at is also reported by Deason et al. (2013) in their proper motion studies of the main-sequence halo stars obtained from the Hubble Space Telescope (HST). This measurement of has broken the mass-anisotropy degeneracy at least out to (K12). To address the problem, practices such as reporting the masses for some arbitrary set of (e.g. Battaglia et al., 2005) or assuming it from simulations (Xue et al., 2008) are a reasonable start. However, to avoid a bias, an approach of marginalizing over all possible values of must be taken. It is worth noting that an independent approach for a mass modeling using halo stars and assuming Jeans Equation is only a good starting point and helpful to make a testable prediction.

One preferred approach (Battaglia et al., 2005; Xue et al., 2008; Kafle et al., 2012) with the Jeans Analysis is to decompose the Galaxy into its dominant components (disk, bulge and halo). Inherent degeneracies among the components is a major concern of this approach and it means assuming a higher disk and/or bulge mass would lower the halo mass and vice-versa. To break the degeneracies, the entire parameter space should be explored. An alternative approach is the tracer mass formalism (Bahcall & Tremaine, 1981; Watkins et al., 2010) which is based on moments. It is a robust technique to estimate the underlying mass of the system provided the density and mass profiles are power laws and the anisotropy is constant with radius, which is certainly not true for the Galaxy.

In this paper, we work towards constructing a holistic model of the Galaxy by combining best available data in hand such as the proper motion of SgrA, the gas terminal velocity in inner of the Galaxy, and kinematics of the large number of halo tracers provided by SDSS/SEGUE. Our aim is to provide the stellar halo number density and kinematic profiles out to the maximum observed distance, the dark matter halo mass and concentration, and the disk and bulge mass parameters. For this we take an approach of fitting parametrized mass models of the Milky Way to observational constraints. This is largely similar to earlier works by e.g. Dehnen & Binney (1998); McMillan (2011); Irrgang et al. (2013) but includes detailed modeling of distant halo stars. Such a model must reproduce known local standard estimates such as the escape velocity (Smith et al., 2007; Piffl et al., 2014), the dark matter density GeV cm (Catena & Ullio, 2010; McMillan, 2011), the total column density (Kuijken & Gilmore, 1989c, a, b, 1991; Bovy & Rix, 2013; Zhang et al., 2013) and the angular velocity of Sun with respect to galactic center (Reid & Brunthaler, 2004; McMillan & Binney, 2010).

We organized the paper as follows; first, in Section 2 we discuss the giants data, outline the selection criteria (for diagnostic see Appendix A) and estimate the distance (for full calculation see Appendix B). In Section 3 we present the halo kinematical profile. In Section 4 we discuss the models for density, anisotropy and potential that are used to fit the kinematics of the halo. and in Section 5 we present our result and discussion. Finally, we summarize in Section 6.

2. Data: Giant Stars

Among the wide varieties of known halo tracers, here we are interested in K giants. These have long been studied (e.g. Ratnatunga et al., 1989; Ratnatunga & Freeman, 1989; Morrison et al., 1990, etc) to probe the distant halo. K giants are brighter, hence, effectively goes deeper. Additionally, they are abundant in number in SEGUE (Yanny et al., 2009), a spectroscopic sub-survey of SDSS. They can therefore supplement the existing catalogs of distant tracers such as the BHB stars (Yanny et al., 2000; Sirko et al., 2004; Xue et al., 2008) and the variable stars (Watkins et al., 2009; Sesar et al., 2011).

2.1. Selection Criteria

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


where and are the extinction corrected magnitudes. The metallicity [Fe/H], and the stellar parameters we use, i.e., surface-gravity and effective temperature are the one labeled in the SDSS DR9 as ‘ADOP’, meaning average of various estimators SSPP 1 pipeline uses, and also are the one 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, we take a slightly conservative cut on to minimize the contamination of dwarfs. Additionally, we apply 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 SSPP provides. Combination of the above given sets of selection criteria yield 5330 candidate giants. For the further diagnostic about selection criteria and candidate giants see Appendix A.

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

2.2. Distance Estimation

Correct distance measurements of the giant stars is critical for studying the kinematics of the halo and also for modelling 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 same as given by Xue et al. (2014). The essence of the exercise is that for each star with some set of observables, say = {, , [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 standard photometric parallax relation to compute the distance for a star. Instead of a single valued number, the set up allows us to compute a full probability distribution or posterior distribution of the distance modulus .

3. The velocity dispersion profile

With the distance modulus posterior distributions and the line-of-sight velocities we now proceed to determine kinematic profiles of the stellar halo, namely the radial , polar and azimuthal velocity dispersion profiles. The line-of-sight velocities we analyze are the one provided by SSPP under the heading “ELODIERVFINAL”. We find that 89% of our sample have uncertainties in radial velocity of . These velocities are in the heliocentric frame of reference and for the purpose of modeling we need to transform them into the one centered at the Galactic center. For this we assume the velocity of the local standard of rest to be an IAU adopted value of ; the motion of the Sun with respect to the local standard of rest to be (Schönrich et al., 2010) and the distance of the Sun from the center of the Galaxy to be .

Had the full space motion of the stars been known, one could measure , and by simply dividing the sample in the radial shells and then computing the second moment of the components of the velocity. Unfortunately, we only know the velocity along one direction, i.e., along the line-of-sight, therefore, obtaining all three dispersion profiles requires a more careful modeling of the halo kinematics.

Given that we are only interested in the dispersion profiles, we consider a gaussian velocity ellipsoid model with rotation about -axis. This model does not require any a priori knowledge of the underlying potential or the tracer density distribution . Generally, the velocity ellipsoid can have a tilt, although it is evident from the recent studies by Smith et al. (2009a) and Bond et al. (2010) that the tilt with respect to the spherical coordinate system for the Galactic halo is consistent with zero so we ignore it. Therefore, we write the velocity distribution as a function of radius as


where the model is given by

the notation represents the gaussian distribution centered at and with dispersion given by


Here, are grid points in radius , which we call nodes. Each of these nodes will have a corresponding value of the velocity dispersion. Thus the final dispersion profile is obtained from a linear interpolation over the nodes.

While the location of the nodes can be fixed arbitrarily, for a more systematic approach, we choose them such that for each node has 500 stars. For , due to fewer stars, we choose them such that each node has 30 stars. This is a non-parametric approach to obtain a kinematic profile and is a useful technique in our case for two reasons. First, we do not have an exact distance but a probability distribution of distance modulus of each star in our sample. Hence, unlike previous studies, e.g., Kafle et al. (2012, 2013) the data cannot be segregated into the radial bins because a star near a bin edge could have some finite probability to be in a neighboring bin. In fact, depending on the distance probability of each star, this approach enables each star to make appropriate contribution to each node where is non zero. Second, in Kafle et al. (2012) undulations were reported in the kinematic profiles of the BHB stars. Hence, we do not want to restrict our analysis by making an assumption about the functional forms of , , and .

In the absence of proper motion information, we marginalize Equation 5 over the tangential velocities, and . The resultant marginalized distribution function (DF) can be expressed as,


Above DF is convolved with distance modulus posterior of each star from Equation B1. The convolution corrects for the spatial selection effect and additionally, also enables to propagate the distance uncertainties to our final estimate of the kinematic profiles.

Finally, to obtain the dispersion profiles we use the likelihood estimation technique based on Markov Chain Monte Carlo (MCMC) random walks. For all experiments below, we use a MCMC algorithm, namely a stretch move as described in Goodman & Weare (2010) to sample from the posterior probability distribution given by our model. Our MCMC walks were run for sufficient autocorrelation time, so as to ensure that the distributions of parameters were stabilized around certain values. Advantage of this method over standard MCMC algorithms such as Metropolis-Hastings is that this method explores the parameter space efficiently and also, produces independent samples with a much shorter autocorrelation time.

The log-likelihood function, , we use is


where sum is over the total number of stars . The MCMC run gives us the posterior distributions of the model parameters at given distances. The values corresponding to the highest likelihood are considered as the best estimates of the model parameters and the uncertainties are computed from the and percentile of the distributions.

Figure 1.— Radial velocity dispersion profile, , of the stellar halo: the black squares with error bars are the values for the giants computed in this paper, the red squares with error bars are the estimates for BHB stars taken from K12, the blue dot with error bar is the measured value for BHB stars from Deason et al. (2012), the magenta dots with error bars are Battaglia et al. (2005) combined estimate for mixed sample of halo objects populating at large and the green dot with horizontal and vertical error bars are the reported value for SEGUE sub-dwarfs by Smith et al. (2009b). The dashed black line is our best fit model for a combined BHB and giants sample (the shaded purple region).

Our profile of the halo giants is shown in Figure 1 by black squares. It starts high at , in the inner region, drops to at distance and then remains flat till . This is consistent with the trend seen in BHB stars (Kafle et al., 2012), shown by red squares. However, note the break in profile of BHB stars is at a radius of , which is slightly smaller than that of the giants. We suspect this to be due to larger distance errors of giants. Distance errors will have the effect of smoothing sharp transitions. At radius , for giants we find that there is a further drop in the reaching as low as at . The magenta data point taken from Battaglia et al. (2005) also shows a low at such a large distance. Similarly, in the range Deason et al. (2012) find a low among their BHB sample (blue data). Similar trends in for different populations reported above is reassuring, since they trace the same gravitational potential. The large error bar in value for the giants sample at is because, close to the Sun and along the galactic pole, the contribution of the radial velocity to the line-of-sight velocity is low.

Figure 2.— Posterior distributions of the velocity dispersions of the giants at kpc: the red, black and blue lines show the probability distribution functions for respectively.

Although we can measure , we are unable to constrain the tangential components of the dispersion, and , for the following reasons. First, with a line-of-sight component of the velocity alone we cannot measure the tangential dispersions at distances . Second, the distance uncertainty of our sample of giants is large, rendering large uncertainties in the tangential dispersions. However, at the first node kpc the tangential velocity contributes significantly to , which makes it possible to compute and , see Figure 2. Here we find that the halo is radial, , which is in agreement with the previous results using the subdwarfs (Smith et al., 2009b), the main-sequence stars (Bond et al., 2010) and the BHB stars (Kafle et al., 2012) at a similar distance.

The issue we ignore in this study is the effect of substructures on kinematics. In K12, it was shown that the effect of the two most dominant structures in the halo, the Sagittarius stellar stream and the Virgo over-density, is negligible.

4. Fitting the model

We now proceed to modeling the profile of the halo in order to probe the Galactic potential, . For this we rearrange the Jeans equation (1) as


We solve this first order differential equation in numerically, forcing a boundary condition as . For some supplied profile for potential , density and anisotropy , we use the numerical solution for to fit the observed obtained in Section 3 (Figure 1). In the following sections we describe parametric forms of , , and that are used in our model, for a quick reference see Table 1.

Physical Quantity Model Comment

Stellar Halo

Density is the break radius, is the truncation radius and is the scale length of fall. Assuming, the logarithmic slope at is continuous gives
Anisotropy Observed values are taken from the literature (Kafle et al., 2012; Deason et al., 2013), 25 is maximum observed distance in kpc, is the velocity anisotropy outside radius .

Dark Halo

NFW potential (Navarro et al., 1996) , and are the virial mass, concentration and the virial radius respectively. is an over-density of the dark matter compared to the average matter density, is the Hubble constant and is the matter density of the universe.


Miyamoto & Nagai (1975) potential , and are the scale length, scale height and mass respectively.


Spheroidal bulge (Binney & Tremaine, 2008, Equation 2.207 a and b) density is a density normalization, which is a function of . is a scale length. Values for remaining parameters , and are assumed to be 1  kpc , 1.8 and 0.6 respectively and are adopted from §2.7 page 111 of Binney & Tremaine (2008).
Table 1Model prescription of the assumed components of the Galaxy.
Galaxy components Parameter Without prior With prior
Stellar Halo ( kpc )
( kpc )
( kpc )
Dark Halo
Disk ()
( kpc )
( kpc )
Bulge ()
( kpc )
Angular velocity ()
Table 2Best-fitting values of the model parameters

4.1. Density

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 () and a sharp fall-off outside. For example analyzing the main-sequence turnoff stars, Bell et al. (2008) conclude that a power-law, , with index of for the inner region and 4 in the outer region with a break at is a reasonable representation. Similar conclusions were also made by Watkins et al. (2009) in their studies of RR Lyrae stars out to and by Deason et al. (2011) for BHB stars out to . Apart from the SDSS survey, the study by Sesar et al. (2011) of the main-sequence turnoff stars obtained from Canada-France-Hawaii Telescope Legacy Survey suggest slightly shallower fall of 3.8 beyond the break radius . More recently, Sesar et al. (2013) using variable stars suggest a broken power law but with much smaller break radius of . In an agreement with Watkins et al. (2009); Deason et al. (2011), here we adopt a double power-law fit to the halo density with inner-slope of and the outer-slope of with break at radius . As we will discuss later, the drop seen in the profile at seems likely to be a consequence of the break in the density. Hence, to further investigate this issue we keep the break radius () as a free parameter of our model.

The outer-most, , 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 comparatively larger number of outer-most halo stars, we investigate as to what does the decline of after tells us about the density profile. For this, we investigate a truncated model of the outer halo with a truncation radius . A similar approximation has also been made by Dehnen et al. (2006), but with a forced sharp truncation at . 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 and the scale length of the fall . The first row in the Table 1 shows the form of the density profile that we adopt. Note, we assume the logarithmic slope at to be continuous. This gives , and the slope at as .

4.2. Anisotropy

The velocity anisotropy () is the most uncertain quantities that enters the Jeans equation, and only recently has been measured directly. The full phase-space studies of the sub dwarfs in Smith et al. (2009b) and main-sequence stars in Brown et al. (2010) suggest that the halo within is radial with . Recently, K12 studied the line-of-sight velocity of BHB stars and found that at similar distance range () the halo is radial . Moreover, K12 also found that the profile of the stellar halo has features which include a tangential dip () at . In the range , 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 the halo is isotropic . For , so far there has been no direct measurement of . Therefore, in this regime we have no choice but to assume some model for . A trend that rises from 0 and attains a positive value at large distances has been reported (see Figure 13 of K12) for the 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 , is constant and equal to value. However, in between the maximum observed radius and the distance , is assumed to be linear. This is more physical than introducing an abrupt transition like a step-function, In other words, determines the slope of the profile joining and . Here , we assume to be . This transition is also consistent with the results from the simulation. However, in both the above mentioned simulations the transitions cease at much closer radius than we assume here. Formula shown in the second row of Table 1 summarizes our model for the anisotropy profile.

4.3. Potential

The next quantity required for the modeling is the potential, . 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 motivated Navarro-Frank-White (NFW, Navarro et al., 1996) profile. The formulas are given in Table 1.

The Galactic dark matter halo assumed to follow an NFW profile, is characterized by two parameters, the virial mass and the concentration ( row of Table 1). Here, the over-density of the dark matter compared to the average matter density, , is considered to be 340 (Bryan & Norman, 1998). The values for Hubble constant and the matter density of the universe are taken from Komatsu et al. (2011) 2; in other words, we assume the mean density is times the cosmological critical density.

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

The Galactic bulge is assumed to be spheroidal. This provides a reasonable axisymmetric approximation for bulge seen from COBE/DIRBE near-infrared data. It is also similar to the axisymmetric approximation of 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 Equation 2.207 a and b in Binney & Tremaine (2008). To compute the force generated by such a spheroidal system we use Equations 2.129 a and b of Binney & Tremaine (2008). The values for some of the bulge parameters are kept fixed, such as, oblateness parameter and power-law index . These were adopted from §10.2.1 of Binney & Merrifield (1998). We found that at distance greater than 4 times the scale length , 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 , which we average over the spherical shells. The average force exerted on a unit mass at given radius due to all the three components of the Galaxy is modeled by

Figure 3.— Effect of density and anisotropy parameters in the model: (a) shows an effect of assumed density power-laws for the case of constant , and and (b) shows an effect of assumed profile and its parameter and for the case of double power-law density and constant and . In all cases bulge and disk parameters are kept fixed as , , , , , , and .

The dominant contributor to the overall potential of the galaxy in the innermost region is the bulge, in 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. So, we keep the parameters free and allow the data to resolve the degeneracies itself.

Before proceeding with the fit, we first study the effect of density and anisotropy profiles in the model. The model is obtained by substituting the above described density, anisotropy and potential profiles in Equation 9. In Figure 3a, we demonstrate the role of the assumed density profile for the case of , and . The solid (dotted) line is the case of a single power-law with slopes 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, break being at . In the figure the 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 . Furthermore, the dashed-dotted line in the figure is same as in the case of dashed line but with an exponential fall-off beyond with parameter . Note the for the dashed-dotted line and dashed line cases are same out to but beyond that the dashed-dotted line declines further due to added exponential fall-off in the density profile. Since the slope in the latter case is a function of the transition is smoother than one near the first break . For , decreases gently. This suggests that the break in the profile is a result of the break, , in the density distribution.

In Figure 3b, we demonstrate the effect of the underlying profile for the case of , , density power law index of 2.4 (inner region), 4.5 (outer region) and . The plotted here are taken from second row of Table 1. The dotted line and dashed-dotted lines show with and and respectively. These two runs can be compared to see the effect of on the model. Clearly, the chosen value of that determines the distance at which saturates and attains constant value, is just the measure of the slope of the transition. Also, as expected, in the case of slope is smaller than in the case of as it attains the at larger r. The effect of seems minimal in the overall profile. Hence, we keep it fixed at for the rest of the analysis. Furthermore, the above two runs can be compared with the dashed line to see the effect of adopted which is 0.5 in the former runs and 0.9 in the latter one. As expected, larger means the shifts up and vice versa. The can systematically shift the 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 (), truncation radius (), truncation softening () and maximum anisotropy (); dark matter halo concentration () and the virial mass (); disk scale lengths ( and ) and mass (); and bulge scale length () and mass (). We consider flat priors for all the parameters in the following range: , , , , , , , , , , and .

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

where is the set of model parameters we explore. The last term is the prior on . The first term is


In this the model is given by Equation 9. The probability is the posterior distribution of parameter at the th 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 at each 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 in the inner most region of the Galaxy , 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 for latitude it is possible to derive a measurement of the rotation curve of the inner Galaxy. Assuming that the Galaxy is axisymmetric, as a function of is given by


where, = . Now we define our second term in likelihood as


There will be effects of non-axisymmetry of the Galaxy and non-circular motion of the ISM on the . To take into account this effects, following Dehnen & Binney (1998), we assume and avoid the region affected by the bar by only using data with . The data with assumed uncertainties are shown with black points in Figure 4

Figure 4.— Terminal velocity as a function of galactic longitude : the points are the data taken from Malhotra (1994, 1995). The error bars of shown are introduced to allow for non-circular motions. The over plotted line is our best fit model resulted from our final MCMC run corresponding to the Figure 6.

An additional prior we impose is on . There is a wide variation in claims about at ranging between (Olling & Merrifield, 1998) to (Méndez et al., 1999). Many of these claims depend on the assumed 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 . As a summary of all these works, we assign a prior with a uniform distribution of


The range in of corresponds to a range in of at .

Figure 5.— Angular velocity () at : the black histogram shows the posterior distribution of obtained from the MCMC run for the case with uniform prior. The best fit value of is . Blue dotted line is the normal distribution with a mean value of and a dispersion of 0.2 obtained for V and from Reid & Brunthaler (2004), which we later use as a prior.
Figure 6.— The joint likelihood and the marginal posterior distributions of the model parameters obtained from the MCMC exploration of the combined sample of the halo giant and BHB stars: the labels along the horizontal and vertical direction tell the name of the parameters, i.e., from the left to the right they are the NFW dark matter halo concentration , virial mass in , density break radius in  kpc , truncation radius in  kpc , truncation softening parameter , anisotropy , disk mass in , disk scale length in  kpc , disk scale height in  kpc , bulge mass in , and bulge scale length in  kpc . Histograms at the top of each column show the posterior distribution of the model parameters named at the bottom of the column whereas the heat maps depict the joint likelihood distribution of two parameters named immediately below and on the left-most end of the same row. Black lines mark the confidence contours. The values in the title of histograms present the best fit estimates.
Figure 7.— Concentration ()–virial mass () contours for : red and green contours are for the giants sample, blue contour is for the BHB sample, and black contour is for a combined sample of BHB and giant stars in a separate distance ranges as labeled in the figure. Closed lines depict region. A white star denotes the best fit estimate and pixel plot shows a 2D posterior distribution corresponding to the black contour. The black dashed line demonstrates a typical relation predicted by CDM dark matter simulation.
Figure 8.— The circular velocity curve of the Galaxy: the dotted, dashed-dotted and dashed lines are the circular velocity curves along the meridional plane, z=0, for the oblate bulge, Miyamoto-Nagai disk and NFW halo respectively. The radius is the distance in the Galactic plane. The individual curves are constructed from the best fit estimates of the model parameters. The solid line shows the resultant circular velocity curves due to all the three components of the Galaxy.
Figure 9.— Escape velocity of the Galaxy: the observed escape velocity of the Galaxy shown as a function of the galactocentric distance . The at is shown by a black dot with error bar. Red stars shows the total velocities for the named Milky Way classical satellite galaxies adopted from Table 1 of Pawlowski & Kroupa (2013).
Figure 10.— Cumulative mass of the Galaxy: the shade shows the observed mass of the Galaxy as a function of the galactocentric distance . The black dot with error bar is Law et al. (2005) estimate of the Galaxy mass within obtained by modeling the Sagittarius dwarf spheroidal galaxy tidal streams whereas the diamond point is mass within obtained by Newberg et al. (2010) by modeling the Orphan Stream.

We have the run for two different tracers, i.e., giant and BHB stars, labelled with red and black in Figure 1. While the run of the giant data is measured out to a larger distance than of the BHB stars , the giant stars have comparatively larger uncertainties in and than BHB stars. Importantly, the 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 : the BHB data from K12 in the range and the giants data in the range . Note, our model for the (see Figure 3) does not predict the flattening-out of the profile in the inner-region . However, the first two values in the figure for the BHB sample show a clear flattening. This is something that needs to be investigated in future, presently we ignore these data points.

During a MCMC run, for every proposed set of values of model parameters, particularly one defining the potentials, we get a prediction for (shown in Figure 5 by a hatched histogram). We are able to constrain the . Interestingly, this is within the uncertainty range of Reid & Brunthaler (2004) result for obtained using the SgrA’s proper motion. Also, our estimate falls in the prescribed range in McMillan & Binney (2010), who use the proper motions of masers. Our measurement therefore can be taken as an independent measurement of the at . Since our constraint of has larger uncertainties compare to the one obtained using SgrA data, from here on unless otherwise mentioned we use as a prior instead of a uniform distribution.

5. Result and Discussion

Figure 6 displays the marginalized 1D and 2D distributions obtained from the MCMC. For a quick referral the table with the best fit estimates is put alongside the figure. It is also summarized separately in Table 2 for both cases, i.e., with and without SgrA proper motion prior. The reported best fit values are the medians of the posterior distributions of the model parameters whereas the uncertainties quoted are computed from and percentile values. The best fit model obtained after substituting the above estimates is shown, against the actual data, by a black dashed line in Figure 1. It can be seen that the best fit represents the data well. A small rise in data at seems to be a local effect and could be due to the presence of some kind of shell like structure at the given distance. For a proper fitting of such outliers we need to know the underlying run of the data.

The Figure 6 nicely demonstrates correlations that exist among different parameters we consider here. For example, one can observe an expected correlation between , also known as the mass-anisotropy degeneracy. Also, one can see a mild correlation between and ; and etc. An anticorrelation is seen between , and and .

The best fit estimates of the model parameters enable us to construct the rotation curve of the Galaxy (shown in Figure 8). The dashed, dotted, and dashed-dotted lines are the as a function of the cylindrical radius for the halo, disk and bulge respectively whereas the solid line is the resultant curve. Substituting the best fit rotation curve in Equation 12, gives us the best fit terminal velocity curve. In Figure 4 it is plotted alongside the terminal velocity data.

5.1. Properties of disk and bulge

In this paper we use a Miyamoto-Nagai 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 which in turn depends upon the chosen value of . Overall the disk mass is around and the bulge mass is around . The structural parameters like and are difficult to measure. The bulge mass has a mild dependence on but other than that and has little effect on other parameters (see Figure 6). We find that, is well constrained by the data. This can be seen by the strong anti-correlation between them and a narrow spread around it. Our inability to measure is due to the following two reasons. First, the priors from terminal velocity data and SgrA we use basically provide and this sensitive only to the sum . Second, in our set up the halo kinematics responds to forces averaged in radial shells, which decreases sensitivity to .

Since, and 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 range and scale heights to be roughly in range . Here we use a Miyamoto-Nagai (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 , we get a scale length of about . Fitting the density in vertical direction, in the range near Sun, we get a scale height of . We adopt a uniform prior for in range , which is within the range expected for exponential disks. The value of that we get is also within the range of expectation. For we adopt a uniform prior in range , which is around the value as suggested by Binney & Tremaine (2008).

5.2. Anisotropy in the outer halo

Our best fit estimate for the anisotropy in outer parts is . It is interesting that we are able to constraint , the reason is as follows. In the inner most region, the terminal velocity data and prior on provide information about the bulge and also to some extent disk parameters. In the region the BHBs anisotropy is already known and the kinematics when put in the Jeans Equation provides estimate for and disk parameters. Now, beyond where is introduced, it is in some sense the only unknown.

5.3. Mass and concentration of dark matter halo

We estimate the mass of the dark matter halo to be and the concentration . Corresponding values for the virial radius and virial velocity are found to be and respectively. It can be seen in Figure 7, that there is a strong anti-correlation between and . The upper bound on is not as well constrained as the lower bound. Simulated virialized halos in 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 7 shows one such relation, , adopted from Macciò et al. (2007). We see that the prediction of the simulation for the dark matter halo in range passes through our measurements. However, note that the predictions are for pure dark matter simulations, and does not include baryonic processes like 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 if the concentration we get is typical or atypical of the Milky Way sized galaxies.

5.4. Do the kinematics of the giant and BHB stars result in a consistent Galactic potential?

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., . For BHB sample we estimate , whereas for giant sample we estimate and . The joint-likelihood distributions for the above two data sets are shown by the blue and green contours in Figure 7, respectively. The distributions for both the samples seem to be in good agreement, except for the fact that it is more puffed for giants than BHB sample. We find that this is mainly due to larger uncertainties in the values for the giants in compare to the BHB sample which we verified by swapping the error bars in BHB and giants. To conclude, the kinematics of the two halo population, namely BHB and giants, are consistent with the fact that they both feel the same Galactic potential, at least within the radius .

5.5. Break in slope of stellar halo density profile

Kinematics of halo stars allow us to constrain the density profile of the halo. We had modelled 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 BHB we measure to be and for giants we measure it to be . A reason for the different for these two halo populations could be the uncertainties in the distances, which is larger for the giants than BHBs. Our estimates for the break radius is slightly smaller than as claimed in Deason et al. (2011) or 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 a power law at much smaller radius of . A smaller break radius also complies with the study of SDSS main sequence turn-off stars in Bell et al. (2008) who conclude that the slope of the density profile at should be shallower in comparison to the radius outside this range. Here, a point worth noting is that our estimate of break radius is linked to the kinematic features whereas all the above values from the literatures are inferred from the studies of the spatial distribution.

5.6. What more do we learn from the tracers extending out to the “edge” of the Galactic halo?

To find an answer we now run the MCMC hammer over the halo giant stars spanning . The red contour in Figure 7 shows the corresponding joint distribution, which is found to almost coincides with the green contours obtained for a giant catalog with  kpc . Similarly, a comparison of the blue contour in the figure, which is for the BHB sample within , with black contour which is for a combined sample of BHB and giant stars in distance ranges , show that they are similar. This suggests that the giant data between does not add much to our knowledge of and , i.e., the potential of the Galaxy. This is mainly due to our adoption of a parametrized form for the distribution of dark mater, namely the NFW profile. The NFW profile predicts that for the falls as . Therefore, data that extends out to about two times the scale radius should be sufficient to constrain the two independent parameters and of the NFW profile. However, if one wants to really compute the density out to virial radius, e.g., by non-parametric schemes then kinematic data till the virial radius would be required. In our case the distant giant stars ( kpc), turn out to be useful to probe the density distribution of the tracer population, i.e. the stellar halo. We find and . 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 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.

5.7. Repercussions of the lighter halo

The number of sub-haloes 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 paradigm predicts larger number of massive subhalos for the Milky Way size 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 against the probability of containing three or less than three subhalos with (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 to compute the mass interior to from the center of the halo at which the mean density is 200 times the critical density. We obtain and corresponding concentration . Figure 5 in Wang et al. (2012) suggests that for the mass equal to our there is probability that the Galaxy should host three or less than three subhalos with . Interestingly, from observations it is known that there are only 3 brightest satellites of the Galaxy namely, Small Magellanic Cloud, Large Magellanic Cloud and Sagittarius dwarf have . This at least suggests that there is no discrepancy between the observed number of luminous satellites with and the number predicted by . Furthermore, Vera-Ciro et al. (2013) also concludes that for Milky Way mass , 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 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 using



Our estimate of as a function of the galactocentric radius 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., at . Except Leo II, which seems to be marginally unbound. We can conclude from the Figure 9 that all the given satellites are bound despite our low estimate for the Galaxy mass.

In Figure 10 we present the cumulative mass of the Galaxy. It is computed using the formula for the centrifugal equilibrium