THE GRAVITATIONAL POTENTIAL NEAR THE SUN FROM SEGUE K-DWARF KINEMATICS

Since the works of Oort (1932, 1960), determining the Galactic gravitational potential from the distribution and kinematics of stars has been one of the most instructive problems in Galactic disk study. Comparing this with the mass density distribution of visible material (stars and gas), one can derive the distribution of dark matter in the Galactic disk. Comparisons of such local dark matter estimations with rotation curve (Weber & de Boer 2010) constrain the shape of the dark matter distribution.

Kuijken & Gilmore (1989a, 1989b, 1989c, hereafter KG89a, KG89b, KG89c) pioneered a practical approach to constrain the Galactic potential near the Sun. They used K-dwarfs to estimate the surface mass density of the total gravitating mass and the mass density of dark matter by solving the "vertical" Jeans equation (i.e., in the one-dimensional equation in the $\hat{\bf z}$-direction), which correlates the space density and velocity of tracer stars. They found the identified (i.e., stellar and gaseous) surface density and the total surface mass density of all gravitating matter within |z| ⩽ 1.1 kpc from the Galactic plane near the Sun to be 48 ± 8 M_☉ pc⁻² and 71 ± 6 M_☉ pc⁻², respectively (KG89a, KG89b, KG89c; Kuijken & Gilmore 1991, hereafter KG91). The uncertainties are mainly caused by the measured errors in distance and velocities of the K-dwarfs and the technique of recovering the vertical force K_z.

After KG91, there have been various other determinations of the disk surface mass density based on different stellar tracers: Flynn & Fuchs (1994, hereafter FF94) derived a surface mass density of 52 ± 13 M_☉ pc⁻² for the known disk matter. With the Hipparcos data, Korchagin et al. (2003) focused Σ_{|z| < 350 pc} = 42 ± 6 M_☉ pc⁻² using red giants; Holmberg & Flynn (2004, hereafter HF04) and Bienaymé et al. (2006) used K-giants and red clump stars to estimate the disk surface density again, and the values were Σ_{|z| < 1.1 kpc} = 74  ±  6 M_☉ pc⁻² and Σ_{|z| < 1.1 kpc} ∼ 57–79 ± 6 M_☉ pc⁻², respectively. Recently, Garbari et al. (2011, 2012, hereafter G11 and G12) revisited the problem, drawing on Jeans equation modeling; they used a Markov Chain Monte Carlo (hereafter MCMC) technique to marginalize over unknown parameters and data from the literature to obtain an estimate of $\Sigma _{|z|<1.1\, {\rm kpc}} = 45.5^{+5.6}_{-5.9}\, M_{\odot }\,{\rm pc^{-2}}$ for the baryonic disk mass and of $\rho _{\rm DM} = 0.025^{+0.014}_{-0.013}\, M_{\odot }\,{\rm pc^{-3}}$. The rather broad range in published values implies that better data and techniques are needed to improve the measurement of the density distribution of the Galactic disk.

In this paper, using the large sample of K-dwarfs observed in SDSS/SEGUE (defined in Section 2; Yanny et al. 2009), we re-determine the disk surface mass density and the mass density of dark matter. In a number of aspects, this work follows KG89b, in particular, in constraining the "vertical force" K_z. However, the paper presents a number of new elements over previous studies in addition to the new large data set. We split the data into abundance-selected sub-samples, which provide distinct probes of the same potential, and we simultaneously fit densities and kinematics using the MCMC approach. We also explore how the results depend on the functional forms for K_z.

This paper is organized as follows. In Section 2, we describe the SDSS/SEGUE data; in particular, we emphasize the spatial selection function, the distance determination, and the sub-samples of similar [Fe/H] and [α/Fe]. In Section 3, we describe how we determine the tracer number density and the vertical velocity dispersion, and the practicalities of solving the Jeans equation. The fitting results are shown and discussed in Sections 4 and 5, with a summary in Section 6.

For our analysis we aim for a sample of kinematics tracers that can be found in disk populations of various ages and metallicities, and which cover a distance of 0.2–2 kpc from the Sun. In the Sloan Digital Sky Survey (SDSS; York et al. 2000), K-dwarfs best satisfy these criteria. The Sloan Extension for Galactic Understanding and Exploration (SEGUE) is a sub-survey of SDSS-II (Abazajian et al. 2009), which operated from 2005 August to 2008 July to probe the formation and evolution of our Galaxy. It obtained ugriz imaging of some 3500 deg² of sky outside the SDSS-I footprint (Yanny et al. 2009, and references therein), with special attention being given to scans of lower Galactic latitudes (|b| < 35°) in order to better probe the disk/halo interface of the Milky Way. Overall, SEGUE obtained some 240,000 medium-resolution (R ∼ 2000) spectra of stars in the Galaxy, selected to explore the nature of stellar populations from 0.3 kpc to 100 kpc (Yanny et al. 2009). The seventh data release (DR 7) was the most recent public data release from SDSS-II, which occurred in 2008 October. SDSS-III, which is currently underway, has already completed the sub-survey SEGUE-II, an extension intended to obtain an additional sample of over 120,000 spectra for distant stars that are likely to be members of the outer-halo population of the Galaxy. Data from SEGUE-II have been distributed as part of the eighth public data release (DR 8; Aihara et al. 2011a). The SEGUE Stellar Parameter Pipeline processes the wavelength- and flux-calibrated spectra generated by the standard SDSS spectroscopic reduction pipeline (Stoughton et al. 2002), obtains equivalent widths and/or line indices for more than 80 atomic or molecular absorption lines, and estimates T_eff, log g, and [Fe/H] through the application of a number of approaches (see Lee et al. 2008a, 2008b; Allende Prieto et al. 2008; Smolinski et al. 2011).

2.1. K-dwarf Selection

KG89b suggested that the tracer stars in the present context should have the following properties: (1) they are phase-mixed and in dynamical equilibrium in the Galactic potential; (2) they are sufficiently common in order to provide statistical precision in the result; (3) they can be found to |z| > 1.0 kpc to ensure that the total surface mass density can be measured; (4) their distances can be well determined. Therefore, considering present SDSS/SEGUE observations, K-dwarfs are ideal stars to use as tracers to measure the total disk mass. Compared to the SEGUE G-dwarfs (e.g., Bovy et al. 2012b), they have the advantage that their minimal distance to be included in SDSS/SEGUE is almost twice as close.

Yanny et al. (2009) list the following color and magnitude cuts for identification and targeting of K-dwarf candidates:

$$\begin{equation*} 14.5 \le \, r_0 \, \le 19.0, \end{equation*}$$

$$\begin{equation*} 0.55 \le \,(g - r)_0 \, \le 0.75. \end{equation*}$$

In order to get reliable estimates for the "vertical" (z) component of each star's velocity, we added two more criteria: the existence of a good proper motion measurement and the existence of spectra in the database,

$$\begin{equation*} {\rm error\; of\; proper\; motion} > 0, \end{equation*}$$

$$\begin{equation*} {\rm S/N} \,> 15. \end{equation*}$$

These criteria returned 10,925 candidates in DR 8. However, DR 8 does not list [α/Fe], which is very important in subsequent analysis for K-dwarf candidates. Besides, proper motions of DR 8 above declination δ > 41° are worse than those of DR 7 because of the astrometric calibration in DR 8 (Aihara et al. 2011b). Thus, we use the [α/Fe] and proper motions from DR 7 for all of these candidates. To eliminate K-giants, only stars whose log g > 4.0 are taken in the present work, which in the end leaves 9157 stars. Figure 1 shows the number density distribution of sample stars in [α/Fe] versus [Fe/H] space.

Zoom In Zoom Out Reset image size

2.2. Abundance-selected Sub-samples

To determine the integral surface mass density and hence measure the Galactic potential near the Sun, using the method described in KG89a, KG89b used K-dwarfs as tracers and assumed that the velocity distribution function of the tracer is only a function of vertical energy: f_z(v_z, z) = f_z(E_z). Adopting the observed vertical tracer density profile ν and a set of parameterized gravitational potential models Φ, which describe the disk surface density of baryonic matter and the volume density of dark matter halo, KG89b used a Jeans equation that includes the σ_Rz term,

$$\begin{equation} -\nu \frac{\partial \Phi }{\partial z} = \frac{\partial }{\partial z}\left(\nu \sigma _{zz}^2\right) + \frac{1}{R}\frac{\partial }{\partial R}\left(R \nu \sigma _{Rz}^2\right), \end{equation} \tag{ 1 }$$

where $\sigma _{ij}^2 = \langle v_iv_j\rangle$ is the velocity dispersion tensor, and an Abel transform

$$\begin{equation} f_z(E_z) = \frac{1}{\pi }\int _{E_z}^{\infty } \frac{-{d}\nu /{d}{\Phi }}{\sqrt{2(\Phi -E_z)}}{d}\Phi \end{equation} \tag{ 2 }$$

to predict f_z(E_z) for each Φ. Convolving it with error distribution, they calculated modeled $f^{\rm mod}_z(v_z,z)$ and compared it with spectroscopic samples to find the best potential parameters. In KG89a's method, the σ_Rz term was also modeled according to the assumption that the axis ratio of the stellar velocity ellipsoid remains constant in spherical–polar coordinates.

In our analysis, we will also adopt the Jeans equation approach and an assumption from KG89b: the Galactic potential is the sum of contributions from baryonic matter and from dark matter. However, we will predict the vertical tracer density profile and the vertical velocity dispersion simultaneously, and we introduce another assumption: the σ_Rz term in the Jeans equation is neglected (see also G11 and G12). To illustrate the validity of the second assumption, we calculate the velocity dispersion tensor directly to estimate the orders of magnitude of the σ_zz and σ_Rz terms in Equation (1) (detailed calculations of velocity components are described in Section 3.2). As in the discussions of the asymmetric drift in Binney & Tremaine (2008; Sections 4.8.2 and 4.9.3), in cylindrical coordinates, assuming that $\sigma _{RR}^2$ and $\sigma _{zz}^2$ both decline with R as exp (− R/R_d), and the velocity ellipsoid points toward the Galactic Center (GC), the σ_Rz term in the Jeans equation can be constrained by

$$\begin{equation} \left| \frac{1}{R}\frac{\partial }{\partial R}\left(R\nu \sigma _{Rz}^2\right) \right| \simeq \frac{2\nu }{R_d}\sigma _{Rz}^2 \lesssim \frac{2\nu z}{R_d} \frac{\sigma _{RR}^2 - \sigma _{zz}^2}{R_{\odot }}, \end{equation} \tag{ 3 }$$

and the σ_zz term is on the order of $\nu \sigma _{zz}^2/z_0$, where z₀ is the disk scale height. We calculate the average value of 〈v_Rv_z〉 to roughly estimate $\sigma _{Rz}^2$, and it is smaller than $\sigma _{zz}^2 = \langle v_z^2\rangle$ by a factor of ∼0.02 for our present sample stars. z₀ ≪ R_d, so the σ_Rz term is smaller than the σ_zz term by at least a factor of $2(\sigma _{Rz}^2/\sigma _{zz}^2)\times (z_0/R_d) \lesssim 0.01$. This satisfies the second assumption. It also implies that we can consider the simple problem of solving the Jeans equation for a one-dimensional slab and then the determine gravitational potential and hence density of matter near the Sun (|R − R_☉| = 1.0 kpc and |z| ≲ 1.5 kpc). Specifically, the "vertical" Jeans equation can be written as (Binney & Tremaine 2008)

$$\begin{equation} \frac{d}{dz} [ \nu _{\star }(z) \sigma _z(z)^2 ] = -\nu _{\star }(z) \frac{d\Phi (z,R_{\odot })}{dz}, \end{equation} \tag{ 4 }$$

where Φ(z, R_☉) is the vertical gravitational potential in the Sun's vicinity, i.e., we solve this equation for R = constant. Here, ν_⋆(z) is the vertical number density of the tracer stellar population and σ_z(z) is the vertical velocity dispersion of those tracers.

With the one-dimensional Poisson equation, which relates the potential to the disk vertical mass density and the first assumption, in cylindrical coordinates we obtain

$$\begin{equation} 4 \pi G \rho _{\rm tot} (z,R_{\odot }) = \frac{d^2}{dz^2 }(\Phi _{\rm disk}(z,R_{\odot })+\Phi _{\rm DM}(z,R_{\odot })), \end{equation} \tag{ 5 }$$

therefore

$$\begin{equation} \rho _{\rm tot} (z,R_{\odot }) = \rho _{\rm disk}(z,R_{\odot }) + \rho _{\rm DM}^{\rm eff}(z,R_{\odot }), \end{equation} \tag{ 6 }$$

where ρ_tot(z, R_☉) is the total mass density of disk components ρ_disk(z, R_☉) and $\rho _{\rm DM}^{\rm eff}(z,R_{\odot })$, an effective dark matter contribution that includes the circular velocity term. In the solar neighborhood, it can be written as (G11)

$$\begin{equation} \rho _{\rm DM}^{\rm eff} (z,R_{\odot })= \rho _{\rm DM}(z,R_{\odot }) - \frac{1}{4 \pi G R}\frac{\partial }{\partial R}V_c^2(R_{\odot }), \end{equation} \tag{ 7 }$$

where V_c(R_☉) is the circular velocity at R_☉. As discussed in G11, if the disk scale height is much smaller than the dark matter halo scale length, the dark matter density can be assumed to be constant over the range of |z|, which is not far from the Galactic midplane, that is, ρ_DM(z, R_☉) ≃ ρ_DM(R_☉). Because in the solar neighborhood the rotation curve is almost flat (Binney & Tremaine 2008; Bovy et al. 2012a), $\rho _{\rm DM}^{\rm eff} \simeq \rho _{\rm DM}(R_{\odot })$. Therefore, it is possible to introduce a third assumption: the dark matter density ρ_DM is taken to be constant within the modeled volume.

Because the first derivative of Φ(z, R_☉) is the gravitational force perpendicular to the Galactic plane K_z(z),

$$\begin{equation} K_z (z,R_{\odot }) \equiv -\frac{d\Phi (z,R_{\odot })}{dz}. \end{equation} \tag{ 8 }$$

Combining Equations (4) and (8), one obtains

$$\begin{equation} \frac{d}{dz} [ \nu _{\star }(z) \sigma _z(z)^2 ] = \nu _{\star }(z) K_z(z,R_{\odot }). \end{equation} \tag{ 9 }$$

To obtain the local mass density distribution from Equation (9), we adopt a parameterized form of $K_z(z,R_\odot \mid \boldsymbol {p})$, where $\boldsymbol {p}$ is set from Equation (6) and it constrains the integral surface mass density of baryonic matter Σ_disk = Σ_⋆ + Σ_gas and its scale height z_h, and the volume density of dark matter ρ_DM. Note that we use Σ_⋆ + Σ_gas rather than ρ_disk, because ρ_disk and z_h are degenerate, i.e., one cannot constrain ρ_disk and z_h separately at the same time.

Then for a given tracer sub-population, we proceed as follows:

• 1.  
  We determine ν_{⋆, obs}(z) and σ_{z, obs}(z) in bins of |z| for each of the three tracer populations directly from the observations.
• 2.  
  We model ν_⋆(z) as a simple exponential (Bovy et al. 2012b) of unknown scale height $\boldsymbol {h}_z$. Here $\boldsymbol {h}_z = [h_{z,1}, h_{z,2}, h_{z,3}]$ correspond to the metal-rich, intermediate-metallicity, and metal-poor sub-samples, respectively.
• 3.  
  We put $\nu _{\star }(z\mid \boldsymbol {h}_z)$ and $K_z(z,R_\odot \mid \boldsymbol {p})$ into Equation (9) and solve for $\sigma _z(z\mid \boldsymbol {p},\boldsymbol {h})$.
• 4.  
  For each [$\boldsymbol {p},\boldsymbol {h}_z$], we compare $\nu _{\star }(z\mid \boldsymbol {h}_z)$ and $\sigma _z(z\mid \boldsymbol {p},\boldsymbol {h}_z)$ with ν_{⋆, obs}(z) and σ_{z, obs}(z) by calculating the likelihood function $\mathcal {L}(\sigma _z,\nu _\star \mid \boldsymbol {p},\boldsymbol {h}_z)$.
• 5.  
  We use a MCMC technique to sample the likelihood of the data, given [$\boldsymbol {p},\boldsymbol {h}_z$] in order to find the best parameters.

In this modeling, we parameterize K_z instead of the gravitational potential. For any K_z function form, the gravitational potential can be derived by integration of K_z. Here we will describe the details step by step.

3.1. Tracer Number Density

3.1.1. Distance Estimates and Coordinate System

To obtain z and v_z of our sample stars in the solar vicinity, we need to first estimate their distance. Given the de-reddened color indices ((g − r)₀ and (g − i)₀), the de-reddened r-band magnitude (in SDSS/DR 8, the extinction corrections in magnitudes are computed following Schlegel et al. 1998), and metallicity, the absolute magnitude is estimated by fitting to fiducial color–magnitude (CM) relations, calibrated through star cluster spectroscopy. The fiducial sequences for (g − r, r)₀ and (g − i, r)₀, based on YREC+MARCS isochrones, as described in An et al. (2009), are adopted. Then it is straightforward to determine the distance $D=10^{(({r_0-M_r})/{5})+1}\,{\rm pc}$.

Such distances were estimated by fitting the fiducials for the two colors discussed above separately, yielding an average difference and standard deviation between these two distance estimates of 0.006 ± 0.077 kpc. We use the average distance and the standard deviation from CM diagrams of (g − r, r)₀ and (g − i, r)₀ as the mean distance and its error, respectively.

With D estimated as above and observed Galactic longitude and latitude (l, b), we get the stars' positions in Galactic cylindrical coordinates. We adopt a Galactocentric cylindrical coordinate system in which the $\hat{\textbf {z}}$-axis is toward the north Galactic pole (b ⋍ 90°), that is,

$$\begin{eqnarray} R&=&\sqrt{D^2\cos ^2(b)-2DR_{\odot }\cos (b)\cos (l)+R^2_{\odot }} \nonumber \\ \theta &=& \arctan \left[\frac{D\cos (b)\cos (l)-R_{\odot }}{D\cos (b)\sin (l)} \right]\\ Z&=& D\sin (b), \nonumber \end{eqnarray} \tag{ 10 }$$

where R_☉ = 8.0 kpc is the adopted distance to the GC (Reid 1993). We aim to obtain the total surface mass density distribution enclosed within a slab of height |z| in the solar neighborhood, as in G12 and Bovy & Tremaine (2012, hereafter BT12).

We will treat the dynamics as a one-dimensional problem (the $\hat{\textbf {z}}$-direction) in the following. This approximation in our analysis can be justified by the limited radial extent of our sample, with a median |R − R_☉| = 0.4 kpc, and by its limited vertical extent, with a median |z| = 0.8 kpc. Therefore, |R − R_☉| and |z| are ≪R_☉ and the neglect of any radial gradients and the velocity ellipsoid tilt should be good approximations.

3.1.2. Selection Function

Any dynamical analysis needs accurate knowledge of the spatial distribution of the kinematic tracer. Therefore, we not only need a well-defined spectroscopic sample of K-dwarfs with known fluxes and distances (or heights above the midplane) estimates, we also need to understand their selection function (Bovy et al. 2012d), which is what we describe here.

In this coordinate system, each sub-sample can be divided into eight bins according to R and ΔR = R_{j + 1} − R_j = 500 pc. In each ΔR range, R is considered as a constant, and each sub-sample is also divided into nine bins in the $\hat{\textbf {z}}$ direction. The width of the z bin, Δz ≡ z_{i + 1} − z_i, is Δz = 100 pc. Here j and i are indices of R and z bins, respectively. In each z–R box,

$$\begin{equation} \nu _{\star,s} (R_j,z_i) = \frac{N_s(R_j,z_i)}{V_{{\rm eff},s}(R_j,z_i)}, \end{equation} \tag{ 11 }$$

where s denotes the different abundance-dependent sub-populations, N_s(R_j, z_i) and V_{eff, s}(R_j, z_i) are respectively the star number and the effective volume, which are corrected by a selection function in a (R_j, z_i) box. Because every targeted star of a given (g − r, r)₀ has a possibility to be in one particular sub-sample, each line of sight is a part of the search volume for every abundance sub-sample. In our effective volume calculation, all lines of sight of all K-dwarfs are included, i.e.,

$$\begin{equation} V_{{\rm eff},s}(R_j,z_i) = \sum _{q=1}^{n_s}V_{{\rm eff},q}, \end{equation} \tag{ 12 }$$

where n_s is the total number of lines of sight of a given sub-sample and V_{eff, q} is the search volume of each line of sight.

It is in the calculation of the effective volume that the selection function enters explicitly. Some aspects of the SEGUE selection function are obvious. The apparent magnitude range brackets the possible distances of K-dwarfs, and the most nearby stars will be preferentially the brightest and coldest stars. At a given r₀ and (g − r)₀, the effective survey volume is smaller for lower metallicity (hence less luminous) stars. Moreover, the SDSS targeting strategy implies that stars at lower latitudes have a smaller probability of ending up in the spectroscopic sample, since a fixed number of targets is observed along each line of sight. In general, we need to derive the selection function, which gives the probability that a star of given M_r, (g − r)₀, and [Fe/H] ends up in the samples as a function of D and (l, b).

Each SEGUE pointing, which corresponds to an angular on-the-sky radius of 149 is observed with two plates: a SEGUE bright plate that targets stars with 14.0 < m_r < 17.8 and a faint plate that targets stars with 17.8 < m_r < 20.1 (Yanny et al. 2009). In each plate, the distribution of the spectroscopic sample in a CM box that satisfies the selection criteria is n_spec(g − r, r), while the distribution of all photometric stars within the same plate that satisfy the same CM cuts is n_photo(g − r, r). Integration of these distributions over all CM cuts results in N_spec and N_photo, leading to a plate weight

$$\begin{equation} W_{{\rm plate},k} \equiv \frac{N_{{\rm spec},k}}{N_{{\rm photo},k}}, \end{equation} \tag{ 13 }$$

where k is the index of the plate. This plate weight should be a part of the selection function. Figure 2 presents the cumulative distributions of spectroscopic and photometric samples using the same selection criteria. It is clear that the selection weight for faint plates has an apparent magnitude dependence because the probability of obtaining a good spectrum is dependent on signal-to-noise ratio but not strongly on color. Therefore, we make an approximation for this selection function, that is, it is a function of r₀.

Zoom In Zoom Out Reset image size

This function is defined for all distance moduli, which is one for distance moduli that correspond to the apparent magnitude inside the plates' magnitude range and zero otherwise. We combine the selection function into the calculation of V_{eff, q}, that is, the average metallicity of one particular sub-sample 〈[Fe/H]〉 is adopted to estimate the possible absolute magnitude, M_q, distance, D_q, and height above the midplane, z_q, of a given (r, g − r). Then,

$$\begin{eqnarray} V_{{\rm eff},q} &=& \frac{\pi }{3}\cdot W_{{\rm plate},q}\cdot W_q\cdot \frac{\sin (\theta)}{\sin (b_{q})}\nonumber \\ &&\cdot \frac{1}{2}[\cot (b_{q} - \theta) - \cot (b_{q} + \theta)] \cdot \left(z_{\rm upper}^3-z_{\rm lower}^3\right),\nonumber \\ \end{eqnarray} \tag{ 14 }$$

where θ = 149 is the diameter of the SEGUE (and SDSS) spectroscopic plate (Yanny et al. 2009) and W_{plate, q} = W_{plate, k} in a given plate. Besides,

$$\begin{equation} W_q = \left\lbrace \begin{array}{@{}ll}1,& {\rm for}\,\,z_{i} \leqslant z_q \leqslant z_{i+1} \\ 0,& {\rm otherwise} \end{array} \right. \end{equation} \tag{ 15 }$$

$$\begin{eqnarray} z_{\rm upper}&=&{\rm Min}(z_{i+1}, D_{r,{\rm max}}\times \sin (b_q)) \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} z_{\rm lower}&=&{\rm Max}(z_{i}, D_{r,{\rm min}}\times \sin (b_q)) \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} D_{r,{\rm min}}&=&10^{\frac{r_{\rm min}-M_q}{5}+1} \,{\rm pc} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} D_{r,{\rm max}}&=&10^{\frac{r_{\rm max}-M_q}{5}+1} \,{\rm pc}, \end{eqnarray} \tag{ 19 }$$

where r_min and r_max are the magnitude limits of each plate. For the SEGUE bright plate, r_min = 14.0 and r_max = 17.8, and for the SEGUE faint plate, r_min = 17.8 and r_max = 20.1 (Yanny et al. 2009). In each z_i bin, the error bar of the tracer number density arises from the star count Poisson variance and it is estimated by means of Monte Carlo bootstrapping (Section 15.6 of Press et al. 2007).

Figure 3 presents the resulting (R, z) map of the tracer number density for each sub-population. From top to bottom, those plots correspond to the metal-rich, intermediate-metallicity, and metal-poor sub-samples; from left to right, plots represent effective star number, effective volume, and natural logarithm of number density in each z-R box. In this figure, one can see that the scale height of the metal-rich sub-sample is shorter, but the scale-length of the same sub-population is longer. This result is similar to the studies of Galactic structure of Bovy et al. (2012d). In each z_i bin, we calculate an average number density over all R and then obtain the vertical profile of the tracer number density.

Zoom In Zoom Out Reset image size

3.2. Vertical Velocity Dispersions

The basic kinematic measurements for each star are its proper motions measured along Galactic longitude latitude $\boldsymbol {\mu }_{l,b}$, and its radial velocity v_rad. From these we calculate the spatial velocity of the star and convert it to the GC by adopting a value of 220 km s⁻¹ for the local standard of rest (V_LSR) and a solar motion of (+10.0, +5.2, +7.2) km s⁻¹ in (U, V, W), which are defined in a right-handed Galactic system with U pointing toward the GC, V in the direction of rotation, and W toward the north Galactic pole (Dehnen & Binney 1998). In the Galactocentric cylindrical coordinate system adopted in the present study, W is the vertical velocity v_z. The error in v_z, $\delta _{v_z}$, is mainly from propagation of observed errors in $\boldsymbol {\mu }_{l,b}$ and v_rad and the error in distance. Figure 4 shows the distribution of $\delta _{v_z}$, which shows that the median error is ∼13 km s⁻¹, and the errors are smaller than 25 km s⁻¹ for 95% of the sample stars. Therefore, it is important that the velocity error of the individual stars is properly folded into the estimate of the velocity dispersion. To estimate σ_z(z, [Fe/H], [α/H]), one therefore cannot simply calculate the standard deviation of the observed v_{z, i} values for each z bin of each sub-population. Instead, we use the maximum likelihood technique described in van de Ven et al. (2006, their Appendix A) to estimate the intrinsic velocity dispersion, corrected for all individual velocity errors.

Zoom In Zoom Out Reset image size

In each z bin, the intrinsic velocity distribution of the stars is assumed as $\mathcal {L}(v_z)$. Each stellar velocity v_{z, i} in this bin is the product of $\mathcal {L}(v_z)$ convolved with a delta function, which is broadened by the observed uncertainties $\delta _{v_{z,i}}$ and integrated over all velocities. For all N stars in the bin, the likelihood can be defined by

$$\begin{equation} \mathscr{L}\,(\bar{v}_z,\sigma _z,\ldots) = \prod _{i}^{N}\int _{-\infty }^{+\infty }\,\mathcal {L}(v_z)\frac{{\rm e}^{-\frac{1}{2}\left(\frac{v_{z,i}-v_z}{\delta _{v_{z,i}}}\right)^2}}{\sqrt{2\pi } \delta _{v_{z,i}}}{d}v_{z}, \end{equation} \tag{ 20 }$$

where $\mathscr{L}\,$ is a function of mean velocity $\bar{v}_z$, velocity dispersion σ_z, and possible higher-order velocity moments. The velocity distribution $\mathcal {L}(v_z)$ can be recovered by maximizing the likelihood $\mathscr{L}$. To describe the velocity distribution, we use a parameterized Gauss–Hermite (GH) series (van der Marel & Franx 1993; Gerhard 1993), because it has the advantage that it only requires the storage of the velocity moments ($\bar{v}_z,\,\sigma _z,\,H_3,\,H_4,\ldots$) instead of the full velocity distribution, and the convolution of Equation (20) can be carried out analytically. For more details on the forms of GH series, we refer the reader to van der Marel & Franx (1993, their Section 2.2 and Appendix A). This makes it feasible to apply the method to a large number of discrete measurements and to estimate the uncertainties on the extracted velocity moments by means of the Monte Carlo bootstrap method (Section 15.6 of Press et al. 2007).

Figure 5 shows an example of a vertical velocity dispersion fit in 600 pc < |z| < 700 pc for the three sub-populations. This figure shows comparisons of different fitting methods. Red solid curves are the results of Gaussian fits without considering observed errors. Green and blue curves represent three and four moments of GH fits using the method described above. For all |z| bins of all sub-populations, the vertical velocity dispersion is overestimated by 2–4 km s⁻¹ when observed errors are not taken into account.

Zoom In Zoom Out Reset image size

For the metal-rich, intermediate-metallicity, and metal-poor sub-samples, the measured vertical velocity dispersion σ_{z, obs}(z) are in the ranges 17–26 km s⁻¹, 25–37 km s⁻¹, and 35–42 km s⁻¹, showing a slight increase with |z|. Moreover, σ_{z, obs}(z) also increases with the decline of metallicity and the increase of [α/Fe] (see Figure 7, filled circles). The most immediately comparable analysis of the vertical kinematics in abundance-selected sub-samples is that of the SEGUE G-dwarfs (Liu & van de Ven 2012; Bovy et al. 2012d). The nearly isothermal dispersions in the mono-abundance bins in Bovy et al. (2012d) that correspond to our metal-rich sample range from 15 to 25 km s⁻¹; those that correspond to our intermediate bin range from 25 to 35 km s⁻¹; those in the metal-poor bin range from 40–50 km s⁻¹. This is in good agreement with our K-dwarf kinematics, except perhaps for the metal-poor bin. However, a quantitative comparison is complicated by the fact that the velocity dispersions expected for populations in broader abundance bins depend on the abundance distribution within that metallicity range (see Bovy et al. 2012d)

3.3. Functional Forms for K_z

Given the observed tracer number density and vertical velocity dispersion, we investigate different parameterized K_z(z) forms to solve Equation (4) and derive the corresponding expression for σ_z(z). K_z(z) has contributions from both the baryonic components and the dark matter. According to the discussion of KG89b, all K_z(z) forms should satisfy the following properties: (1) the halo and disk contributions are degenerate when |z| ≪ scale height; (2) the disk contribution can be approximated by a thin mass sheet when |z| ≫ scale height. We considered five different parameterized forms of K_z(z), two of which we discuss below, while the other three are given in the Appendix.

• 1.  
  KG89 model

  $$\begin{equation} K_z = 2 \pi G \left[\Sigma _{\star } \frac{z}{\sqrt{z^2+z_h^2}}+\Sigma _{\rm gas} \right]+ 4 \pi G \rho _{\rm DM} z. \end{equation} \tag{ 21 }$$
• 2.  
  Exponential model

  $$\begin{eqnarray} &&K_z = 2 \pi G \left\lbrace \Sigma _{\star } \left[1-\exp (-\frac{z}{z_h})\right]+\Sigma _{\rm gas}\right\rbrace + 4 \pi G \rho _{\rm DM} z.\nonumber \\ \end{eqnarray} \tag{ 22 }$$

Here, Σ_gas = 13.2 M_☉ pc⁻² is the surface mass density of gaseous components which is taken from the disk mass model of Flynn et al. (2006).

As we stated before, the trace number density profile of each sub-population is modeled by a single exponential,

$$\begin{equation} \nu _{\star,s} (z) = \nu _{0,s}\times \exp \left(-\frac{z}{h_{z,s}}\right), \end{equation} \tag{ 23 }$$

where s is the index of the sub-population. For the exponential model, the corresponding vertical dispersion profile of each sub-population can be derived analytically:

$$\begin{eqnarray} \sigma _{z,s}^2 (z)&=& 2 \pi G h_{z,s}\left\lbrace \Sigma _{\star }\left[1-\frac{z_h}{h_{z,s}+z_h}{\rm exp}\left(-\frac{z}{z_h}\right)\right]+ \Sigma _{\rm gas}\right\rbrace\nonumber \\ && + 4 \pi G \rho _{\rm DM} h_{z,s}(z+h_{z,s}). \end{eqnarray} \tag{ 24 }$$

For the other K_z forms we infer σ_{z, s} through numerical integration.

Then it is possible to fit σ_{z, s}(z) and ν_{⋆, s}(z) simultaneously and obtain the parameters described below.

3.4. Obtaining the PDFs for the Model Parameters

The predictions for $\sigma _{z,s}(z\mid \boldsymbol {p},\boldsymbol {h}_z)$ are based on six parameters (i.e., $\boldsymbol {p}$=[Σ_⋆, z_h, ρ_DM] and $\boldsymbol {h}_z$ = [h_{z, 1}, h_{z, 2}, h_{z, 3}]) for each K_z(z, R_☉) model and the vertical tracer populations, which may have considerable covariance. At the same time, we have a parameterized model for the vertical tracer density profile $\nu _{\star }(z\mid \boldsymbol {h}_z)$, and we want to derive constraints on the gravitational force, marginalized over $\boldsymbol {h}_z$. To explore the parameter space efficiently and to account for parameter degeneracies, we use an MCMC approach to sample the likelihood function, which provides a straightforward way of estimating the probability distribution functions (PDFs) for all the model parameters and their degree of uncertainty (Section 15.8 of Press et al. 2007).

The estimates σ_z(z) and ν_⋆(z) in each |z| bin are independent, therefore the logarithm of the likelihood of the data given the parameters $[\boldsymbol {p},\,\boldsymbol {h}_z]$ can be written as

$$\begin{eqnarray} \ln \mathcal {L}&= & -\sum _{m=1}^{M}\ln (2\pi \epsilon _{\sigma _{z,{\rm obs}}}\, \epsilon _{\nu _{\star,\rm obs}})_m\,- \sum _{m=1}^{M}\frac{1}{2} \nonumber \\ &&\times \left[\frac{\sigma _z(\boldsymbol {p},\boldsymbol {h}_z)-\sigma _{z,{\rm obs}}}{\epsilon _{\sigma _{z,{\rm obs}}}}\right]^2_m - \sum _{m=1}^{M}\frac{1}{2}\left[\frac{\nu _{\star }(\boldsymbol {h}_z)-\nu _{\star,{\rm obs}}}{\epsilon _{\nu _{\star,\rm obs}}}\right]^2_m, \nonumber \\ \end{eqnarray} \tag{ 25 }$$

where M is the total number of |z| bins, and $\epsilon _{\sigma _{z,{\rm obs}}}$ and $\epsilon _{\nu _{\star,\rm obs}}$ are the errors of observed σ_{z, obs} and ν_{⋆, obs}, respectively.

In order to reduce the number of parameters by eliminating ν_{0, s} from Equation (23), we calculate ν_{0, s} for each sub-population by

$$\begin{equation} \frac{\partial \chi ^2(\nu)}{\partial \nu _{0,s}} = 0, \end{equation} \tag{ 26 }$$

yielding

$$\begin{eqnarray} &&\nu _{0,s} = \left. \sum _{m=1}^M \frac{\exp (-z_m)}{\epsilon _{\nu _{\rm obs},s,m}}\nu _{{\rm obs},s,m} / \sum _{m=1}^{M}\frac{\exp (-z_m)}{\epsilon _{\nu _{\rm obs},s,m}}\exp \left(-\frac{z_m}{h_{z,s}}\right) \right..\nonumber\\ \end{eqnarray} \tag{ 27 }$$

All other parameters are estimated by running a MCMC chain of typically 50,000 steps, after a burn-in period. In the Markov chain, the parameters of the gravitational model and scale heights of tracer densities [$\boldsymbol {p}_n$, $\boldsymbol {h}_{z,n}$] are chosen from a multivariate normal distribution of [$\boldsymbol {p}_{n-1}$, $\boldsymbol {h}_{z,n-1}$], with the covariance of [$\boldsymbol {p}$, $\boldsymbol {h}_z$] chosen to yield a total acceptance rate in the MCMC chain of about 20%–30%. We use the results of all steps to represent [$\boldsymbol {p}$, $\boldsymbol {h}_z$] and their distributions by counting the rejected steps as repeat ones. After sampling the PDF for parameters for each K_z model, $\chi ^2_{\rm tot}$ between the observed values and model predicted ones (σ_{z, mod} and ν_{⋆, mod}) is calculated to explore how the models from the different K_z families compare:

$$\begin{eqnarray} && \chi ^2_{\rm tot} = \chi ^2_{\sigma _z} + \chi ^2_{\nu _{\star }} = \sum _{s=1}^3 \nonumber \\ &&\quad \times \left[\sum _{m=1}^M\left(\frac{\sigma _{z,s,{\rm mod}} - \sigma _{z,s,{\rm obs}}}{\epsilon _{\sigma _{z,s,{\rm obs}}}}\right)^2_m+\sum _{m=1}^M\left(\frac{\nu _{s,{\rm mod}}-\nu _{s,{\rm obs}}}{\epsilon _{\nu _{s,{\rm obs}}}}\right)^2_m\right]. \nonumber \\ \end{eqnarray} \tag{ 28 }$$

The above procedure results in PDFs for the K_z and ν_⋆ parameters, which are presented in Figures 6–11. The best-fitting results and recovered parameters for all K_z gravitational force law models are listed in Table 1.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Vertical Tracer Density Profiles

4.2. Vertical Dispersion Profiles

4.3. The Vertical Gravitational Force K_z

Figure 8 shows K_z(z) implied by our fits using the KG89b and exponential families for the "vertical force law." The long-dashed lines show K_z(z) for the most likely parameters, and the band of gray lines shows a 1σ sampling of the PDF for K_z(z). Both the KG89b and the exponential families of K_z(z) show a number of generic features in Figure 8: they start out with a finite value for small |z|, as we have fixed a prior contribution from a thin layer of the cold gas with ∼13 M_☉ pc⁻²; then K_z rises steeply to ∼300 pc (reflecting the mass scale height of the stellar disk) and then flattens out; beyond ∼500 pc the slower rise in K_z reflects the dark matter halo term.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

On the basis of K_z(z), we can derive the total surface density at |z| = 1.0 kpc. We get $\Sigma _{{\rm tot},|z|<1.0\, {\rm kpc}}^{\rm KG89} = 67 \pm 6\, M_{\odot }\,{\rm pc^{-2}}$ for the KG89 model and $\Sigma _{{\rm tot},|z|<1.0\, {\rm kpc}}^{\rm Exp} = 66 \pm 8\, M_{\odot }\,{\rm pc^{-2}}$ for the exponential model. It is clear that Σ_{tot, |z| < 1.0 kpc} is robustly constrained, irrespective of the K_z model. Our modeling also constrains the mass scale height: $245^{+188}_{-245}$ pc for KG89 and $200^{+100}_{-200}$ for the exponential model (see Figure 9). Within these uncertainties, this is consistent with the 180 pc from Hill et al. (1979) by analyzing A- and F-dwarfs and $390^{+330}_{-120}$ pc of Siebert et al. (2003), found by using high-resolution spectral data of red clump stars. Taken together this firms up the picture that there is a dominant mass layer near the disk midplane (presumably baryonic) that is rather flat.

For comparison with literature studies, we also explicitly estimate the total surface density of baryonic matter within |z| < 1.1 kpc, Σ_{baryonic, |z| < 1.1 kpc} = Σ_⋆ + Σ_gas, which we find to be $\Sigma _{{\rm baryonic},|z|<1.1\, {\rm kpc}}^{\rm KG89} = 55 \pm 5\, M_{\odot }\,{\rm pc^{-2}}$ and $\Sigma _{{\rm baryonic},|z|<1.1\, {\rm kpc}}^{\rm Exp} = 54 \pm 8\, M_{\odot }\,{\rm pc^{-2}}$, respectively. These values are slightly higher than 48 ± 8 M_☉ pc⁻² derived by KG89b, but are in agreement with the value of 52 ± 13 M_☉ pc⁻² of FF94 and 53 ± 6 M_☉ pc⁻² of HF04.

Table 1. Fitting Results of the Two Models, with Σ_gas = 13 M_☉ pc⁻² (Flynn et al. 2006)

	KG89 Model (Equation (21))	Exponential Model (Equation (22))
Σ⋆ (M☉ pc−2)	42 ± 6	41 ± 5
ρDM (M☉ pc−3)	0.0064 ± 0.0023	0.0060 ± 0.0020
zh (pc)	$245^{+188}_{-245}$	$200^{+100}_{-200}$
hz, 1a (pc)	259 ± 12	260 ± 15
hz, 2a (pc)	450 ± 26	465 ± 33
hz, 3a (pc)	852 ± 30	910 ± 71
σ0, 1a (km s−1)	15.4 ± 1.3	15.8 ± 1.3
σ0, 2a (km s−1)	23.0 ± 2.0	23.6 ± 2.0
σ0, 3a (km s−1)	34.2 ± 2.2	35.8 ± 2.2
$\chi ^2_{\sigma _z}$	25.23	26.36
$\chi ^2_{\nu _{\star }}$	15.80	14.54
$\chi ^2_{\rm tot}$	41.03	40.90

Note. ^aNumbers 1, 2, and 3 represent the metal-rich, intermediate-metallicity, and metal-poor sub-samples, respectively.

Download table as:  ASCIITypeset image

The local mass densities of dark matter recovered by our models are $\rho _{\rm DM}^{\rm KG89} = 0.0065\pm 0.0023\, M_{\odot }\,{\rm pc^{-3}}$ (0.25 ± 0.09 GeV cm⁻³)⁶ and $\rho _{\rm DM}^{\rm Exp} = 0.0060\pm 0.0020\, M_{\odot }\,{\rm pc^{-3}}$ (0.23 ± 0.08 GeV cm⁻³), marginalized over the other parameters. We explored the parameter degeneracies by MCMC (see Figure 11), which recover the known degeneracy between the surface mass density of stars and the volume mass density of dark matter: lower Σ_⋆ correlates with higher ρ_DM at an approximately constant total surface mass density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have presented an analysis of the vertical Galactic potential at the solar radius, drawing on ∼9000 K-dwarfs from SDSS/SEGUE. In many ways, the analysis followed the Jeans equation approach initially laid out by KG89 and implemented by several other groups since (Siebert et al. 2003, HF04, G12, BT12). In comparison to most previous studies, our analysis has a number of new elements: we have a substantially larger sample than previous analyses; we have taken explicit account of the abundance-dependent spatial selection function of our sample; we have simultaneously fit for several abundance-selected nearly isothermal sub-populations that "feel" the same gravitational potential, and we have matched the kinematics and the spatial distribution simultaneously; we have explored to what extent the choice of the functional form for K_z affects the results.

As laid out in the previous section, we found good constraints on Σ_{tot, |z| < 1.0 kpc} and some constrains on the thickness of the disk mass layer, z_h: Σ_{tot, |z| < 1.0 kpc} = 67 ± 6 M_☉ pc⁻² and z_h ≲ 300 pc, irrespective of the functional form we assume for K_z. This is also among the first studies to find significant (>2σ) constraints on the local dark matter density, ρ_DM = 0.0065 ± 0.0023 M_☉ pc⁻³ (0.25 ± 0.09 GeV cm⁻³). Figure 11 shows the expected degeneracy between the total surface density of visible matter and the dark matter density: a lower Σ_⋆ corresponds to a higher ρ_DM for a given total surface density. However, these uncertainties do not include the systematic errors from misestimated photometric distances. For example, unrecognized binary contamination would make the inferred distances 10% larger, and the distance estimation of An et al. (2009) that we adopted in the present work is ∼9% smaller than that of Ivezić et al. (2008). We explore distance systematics in particular for their implications on the dark matter density: if we systematically change all distances in the input data catalog by ±10% and then go through the analysis steps considered in Section 3, the inferred change of ρ_DM is only ±10%.

Overall, our results on K_z are consistent with earlier studies, as shown in Figure 12, which compares the inferred total surface density Σ_tot(z) with previous studies. The value derived by KG91 is Σ_{tot, |z| < 1.1 kpc} = 71 ± 6 M_☉ pc⁻² and ρ_DM = 0.01 ± 0.005 M_☉ pc⁻³, while the calculation of HF04 was Σ_{tot, |z| < 1.1 kpc} = 74 ± 6 M_☉ pc⁻² and ρ_DM = 0.007 M_☉ pc⁻³. The latter results are in a good agreement with ours.

Zoom In Zoom Out Reset image size

Figure 13 compares our inference for ρ_DM with two recent studies. While it agrees very well with BT12, our value is considerably lower than the ρ_DM inferred by G12, who found a value of $\rho _{\rm DM} = 0.025^{+0.014}_{-0.013}\,M_\odot \,{\rm pc}^{-2}$ ($0.85^{+0.57}_{-0.49}\ {\rm GeV\,cm^{-3}}$) but the two sets of works are statistically consistent at 90% confidence. G12 used a five times smaller sample over a comparable value; the sample had a simple selection function, but was not separated in abundance sub-bins. This may be the main cause of the difference in the dynamic mass density of dark matter. In principle, the Jeans equation is linear, which means that the dynamical inferences should not change whether the tracer population is mixing or not. However, splitting stellar populations according to their abundances can provide stronger constraints on scale heights of the vertical profile of tracers and hence place a stronger constraint on the dark matter density. We use an exponential K_z model (Equation (22)) as an example to make a test and find that the errors on ρ_DM are larger.

Zoom In Zoom Out Reset image size

On the other hand, our estimates are in very good agreement (Figure 13) with the recent determination of ρ_DM by BT12, who re-analyzed the data by Moni Bidin et al. (2012). They used the assumption of ∂V_c/∂R = 0 to correct that model and derived a dark matter density of ρ_DM = 0.008 ± 0.003 M_☉ pc⁻³ (0.30 ± 0.11 GeV cm⁻³).

Figure 13 puts the recent local ρ_DM estimate from K_z constraints into context. The histogram with the gray shading in its center represents the results from G12, the one with black is from BT12, and the one with red is from the present work (hereafter Z12). The BT12 histogram reflects their ρ_DM PDF including systematic errors; for Z12 several of these error sources have been addressed systematically, and the distance uncertainties have been incorporated as a systematic error term in Figure 13. Therefore, the results shown in Figure 13 from Z12 and BT12 should be on a comparable footing. The joint PDF emerging from these three recent experiments is shown as the thick black histogram in Figure 13, indicating ρ_DM = 0.0075 ± 0.0021 M_☉ pc⁻³ (0.28 ± 0.08 GeV cm⁻³).

One can put these "local" estimates of the dark matter density into the context of the expectations from the global fits of the Galactic dark matter halo: with a spherical Navarro–Frenk–White cold dark matter density profile (Navarro et al. 1996) and the parameters of Xue et al. (2008; virial mass $M_{\rm vir} = 0.91^{+0.27}_{-0.18} \times 10^{12}\, M_{\odot }$, virial radius $r_{\rm vir} = 267^{+24}_{-19}\ {\rm kpc}$, and c = 12.0), one would expect that the dark matter density at R_☉ is ρ_DM(R_☉) = 0.0063 M_☉ pc⁻³, which is in very good agreement with our present work.

Taken together with the existing work, our results continue to point toward a picture that the local Galactic potential is dominated by a fairly thin layer of stars and gas plus a dark matter density that is in accord with global fits to the Galactic halo. Specifically, we find no evidence for any significant amount of disk dark matter: our combined Z12+BT12 estimate rules out the best-fit value of G12 of ρ_DM = 0.025 M_☉ pc⁻³ at 8σ and our dynamical measurement of the surface density of baryonic matter of Σ_baryonic = 55 ± 6 M_☉ pc⁻² is in good agreement with direct estimates through star counts and mapping of the local interstellar medium. But note that due to the error bar of G12, our results are formally consistent with theirs; therefore, a disk-like dark matter distribution at a normalization lower than that advocated by G12 cannot be ruled out, and should be tested further, as it contributes a qualitatively different scenario from an only round halo.

It may seem surprising that our derived limits are not much tighter than those obtained by KG89b and subsequent work, based on smaller samples. This is in good part due to the fact that we used far fewer prior constraints on the models. For example, we fit for the mass scale height of the disk (z_h) rather than assume it; similarly, we fit for the scale height of the tracers simultaneously with the kinematic fit; we fit for dark matter density rather than assume it as a prior; and we have explored a range of functional forms for K_z. Marginalizing over these factors apparently leads to a similar error on the parameters than other more restricted analysis with smaller samples.

We have analyzed the K-dwarfs from SDSS/SEGUE to determine the total surface mass density and dark matter density in the local Galactic disk. First, we divide our sample into three sub-populations through their [Fe/H] and [α/Fe]. After considering the spatial selection function, the Galactic vertical number density profiles for different sub-populations are inferred from star counting. Second, we use maximum likelihood and GH series to calculate the vertical velocity dispersion profile of each sub-population. Then different parameterized K_z forms are used to solve the "vertical" Jeans equation. We fit the observed vertical number density and vertical velocity dispersion profiles of the three sub-populations simultaneously, using a MCMC technique to recover the PDFs of the parameters of the "vertical force law." In our results, for each sub-population, the vertical number density is approximately a single exponential profile and the vertical velocity dispersion is nearly isothermal. The scale height of the number density profile and the vertical velocity dispersion increase with the decline of metallicity. Assuming that there is a thin gas layer with Σ_gas = 13 M_☉ pc⁻², we derive a total surface mass density of 67 ± 6 M_☉ pc⁻² at 1.0 kpc from the midplane of which the contribution of all stars is 42 ± 5 M_☉ pc⁻² and we infer a local dark matter density of ρ_DM = 0.0065 ± 0.0023 M_☉/pc⁻³ (0.25 ± 0.09 GeV cm⁻³).

We thank the anonymous referee for helpful comments, Constance Rockosi for the assistance of distance calculation, and Justin Read for discussions and suggestions. L.Z. acknowledges support of NSFC grants 10903012 and 11103034 and from the MPG-CAS student program. H.W.R., G.v.d.V., and J.B. acknowledge partial support from Sonderforschungsbereich SFB 881 "The Milky Way System" (subproject A3 and A7) funded by the German Research Foundation. J.B. was supported by NASA through Hubble Fellowship grant HST-HF-51285.01 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III Web site is http://www.sdss3.org/.

SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, University of Cambridge, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University.

Besides the two models for K_z described in Section 3.3, we explore the other three models to determine the total surface mass density and the mass density of dark matter, with the purpose of checking whether the astrophysical inferences depend on the choice of these functional forms:

• 3.  
  Error function model

  $$\begin{equation} K_z = 2\,\pi \,G\left[\Sigma _{\star }\,{\rm erf}\left(\frac{z}{z_h}\right) +\Sigma _{\rm gas} \right]+ \,4\,\pi \,G\,\rho _{\rm DM}\,z. \end{equation} \tag{ A1 }$$
• 4.  
  General model I. This form is an extension of the KG model. We assume K_z as of the following form:

  $$\begin{equation} sK_z = 2\,\pi \,G\left[\Sigma _{\star }\left(\frac{z^\beta }{z^\gamma +z_h^\beta }\right)^ {\frac{1}{\beta }}+\Sigma _{\rm gas}\right]+\,4\,\pi \,G\,\rho _{\rm DM}\,z. \end{equation} \tag{ A2 }$$
• 5.  
  General model II. To reduce the number of free parameters, we set β = γ and get:

  $$\begin{equation} K_z = 2\,\pi \,G\left[\Sigma _{\star }\,\left(\frac{z^\beta }{z^\beta +z_h^\beta }\right)^ {\frac{1}{\beta }}+\Sigma _{\rm gas}\right]+\,4\,\pi \,G\,\rho _{\rm DM}\,z. \end{equation} \tag{ A3 }$$

We followed the same parameter estimation approach for these models and calculated the total χ² between observed σ_{z, obs} and ν_⋆ and the model predictions. The error function model and the general model II yield similar results for the surface density of visible matter and dark matter density as the KG89 and exponential models discussed in the main text. Only the general model I, where the exponents γ and β are independent, yields a degree of parameter degeneracy that makes the resulting K_z(z) difficult to interpret. But overall this confirms that, among physically plausible families of K_z, the particular choice matters little.