ON THE LOCAL DARK MATTER DENSITY

The observed flatness of the Milky Way's circular-velocity curve at Galactocentric distances larger than 20 kpc (e.g., Xue et al. 2008) shows that the visible Galactic disk is embedded in a massive dark halo. The disk is composed of gas and stars (baryons), while the dark halo is believed to be composed of non-baryonic matter of unknown nature. Despite the dominance of the dark halo in the outer parts of the Milky Way, it remains unclear from direct measurements whether there is any need for a substantial amount of dark matter to explain the circular-velocity curve interior to the solar radius, R₀ = 8 kpc. Quantitatively, the fraction of the total radial force at the solar radius that is due to the disk could be as high as 90% (Sackett 1997). A promising avenue for constraining the local density of dark matter is through a determination of the dependence of the gravitational potential on distance above the midplane of the disk ("height"), from measuring the kinematics of stars (e.g., Kapteyn 1922; Oort 1932; Bahcall 1984). However, a major obstacle is that the uncertainty in the amount of baryonic matter in the disk makes it hard to determine the relative contributions from dark and baryonic matter to the density near the midplane. Studies limited to heights ≲ 150 pc are consistent with no dark matter near the Sun, but they cannot exclude the amount of dark matter expected for a standard dark matter halo (≈0.01 M_☉ pc⁻³ = 0.38 GeV  cm⁻³; Crézé et al. 1998; Holmberg & Flynn 2000).

The contributions from baryonic and dark matter can be disentangled by measuring the gravitational potential out to larger heights. At heights of several times the disk thickness, the dark halo and the baryonic disk contributions to the potential have a different vertical dependence (e.g., Kuijken & Gilmore 1989; Garbari et al. 2011). In particular, most of the disk mass lies below ∼1.5 kpc, so above this height the disk contribution to the integrated surface density Σ(Z) ≡ ∫^Z_{− Z}dz ρ(z) is roughly constant with height and the disk potential varies as Φ∝|Z|, while the dark halo contributes Σ(Z)∝|Z| and Φ(Z)∝Z². Thus, any measured increase in the surface density at |Z| ≳ 2 kpc must be due to the dark halo, and the expected increase is ∼20 M_☉ pc⁻² for each kpc (for the standard value of the local dark matter density of 0.01 M_☉ pc⁻³; see above). Determinations of the surface density at ∼1 kpc from the plane typically find values of Σ(1 kpc) = 70–80 M_☉ pc⁻² (Kuijken & Gilmore 1991; Siebert et al. 2003; Holmberg & Flynn 2004), while the baryonic contribution is estimated to be around 50–60 M_☉ pc⁻² (Holmberg & Flynn 2004; Bovy et al. 2012a). Thus, these studies are consistent with the expected dark matter density, although they do not go to large enough heights to detect the dark matter unambiguously.

A recent study of the kinematics of stars with heights 1.5 kpc ≲ |Z| ≲ 4 kpc by Moni Bidin et al. (2012b, MB12 hereafter) claims to have measured the surface density at these heights in a model-independent manner and found it to be constant, such that the local density of dark matter must be smaller than 10⁻³ M_☉ pc⁻³ (<0.04 GeV  cm⁻³). If true, this would imply that dark matter is less abundant in the solar neighborhood than expected, by at least an order of magnitude. This would also shift experimental limits on the cross section for elastic scattering between dark matter particles and baryons by an order of magnitude, although this is far less than the uncertainty in the predicted cross section.

In this paper, we show that the analysis used by MB12 is flawed. The main error is that they assume that the mean azimuthal (or rotational) velocity $\bar{V}$ of their tracer population is independent of Galactocentric cylindrical radius R at all heights (i.e., $\bar{V}(R,Z)=\bar{V}(Z)$). This assumption is not supported by the data, which instead imply only that the circular speed V_c is independent of radius in the midplane (e.g., Gunn et al. 1979; Feast & Whitelock 1997). In the solar neighborhood, the circular speed is larger by ≳ 35  km s⁻¹ than the mean azimuthal velocity for the warm tracer population used by MB12, a phenomenon known as asymmetric drift. The asymmetric drift is expected to vary with R—although this variation cannot be measured for the sample of MB12 as the data do not span a large enough range in R—so the assumptions that $\bar{V}$ and V_c are independent of radius are not compatible. In the absence of a measurement of the (R, Z) dependence of $\bar{V}$ for the tracers, we show that the assumption of an R-independent $\bar{V}$ at all heights Z is highly implausible. By using instead the data-driven assumption that the circular speed is independent of radius at Z = 0, we demonstrate that the measurements and analysis of MB12 are fully consistent with the standard estimate of the local dark matter density, approximately 0.01 M_☉ pc⁻³, and indeed provide the best available direct measurement of this quantity.

The outline of this paper is as follows. In Section 2, we discuss the assumption made by MB12 that the mean azimuthal velocity of the tracer population is independent of R. In Section 3, we derive the surface density as a function of Z using the assumption that the circular-velocity curve is flat in the midplane and discuss the effect of relaxing this approximation. In Section 4, we calculate the surface density using the data of MB12 but using the approximation of constant circular speed and show that this leads to a surface density at 1.5 kpc < |Z| < 4 kpc that increases with height at the rate expected for a dark halo. Our conclusions are in Section 5.

We follow MB12 in defining U, V, and W as the radial, azimuthal, and vertical velocities in cylindrical coordinates and the inertial Galactocentric reference frame. Note that this definition is non-standard, since (1) U, V, and W are normally defined with respect to the local standard of rest; (2) U is typically the velocity toward the Galactic center; and (3) our conventions require either that the coordinate system is left-handed or that the positive z-axis points to the south Galactic pole. However, for the purpose of this paper these distinctions do not matter. We use R₀ = 8 kpc, as do MB12, and local circular speed V_c = 220  km s⁻¹ throughout this paper, assumptions that are consistent with the latest measurements (Bovy et al. 2009). We abbreviate the surface density at the solar radius as Σ(Z) ≡ Σ(R₀, Z) = ∫^Z_{− Z}dz ρ(R₀, z).

The analysis of MB12 uses the Poisson equation to calculate the surface density as a function of height Z at the solar radius R₀ using the radial force F_R and the vertical force F_Z (Equation (18)). To estimate F_R, MB12 make use of the radial Jeans equation

$$\begin{eqnarray} F_R(R,Z) &=&\, {-}\frac{\partial \Phi (R,Z)}{\partial R} = \frac{1}{\nu }\,\frac{\partial \left(\nu \sigma _U^2\right)}{\partial R} +\frac{1}{\nu }\,\frac{\partial \left(\nu \sigma ^2_{UW}\right)}{\partial Z}\nonumber\\ && + \frac{\sigma _U^2 - \sigma _V^2 - \bar{V}^2}{R}, \end{eqnarray} \tag{ 1 }$$

where Φ is the gravitational potential, ν is the tracer-density profile, σ²_U and σ²_V are the radial and azimuthal velocity dispersions squared, σ²_UW is the off-diagonal radial–vertical entry of the dispersion-squared matrix, and $\bar{V}$ is the mean azimuthal velocity; all of these quantities are functions of R and Z. We have assumed a steady state such that time derivatives vanish and the mean radial motion $\bar{U}$ is zero. In simplifying Equation (1), MB12 use their assumption VIII: "The rotation curve is locally flat in the volume under study," which MB12 express as

$$\begin{equation} \frac{\partial \bar{V}(R_0,Z)}{\partial R} = 0. \vspace*{6pt} \end{equation} \tag{ 2 }$$

Moni Bidin et al. (2010) made the same assumption, although this was not explicitly stated. We show in this section that this is an unreasonable assumption.

The mean azimuthal velocity of a population of stars differs from the circular velocity due to the asymmetric drift. This offset arises because both the density of stars and the velocity dispersion typically decline with radius. This means that more stars with guiding centers at R < R₀ are passing through the solar neighborhood than stars with guiding centers R > R₀; the former are on the outer parts of their orbits, where their azimuthal velocity is less than the circular velocity. Equation (1) is typically used to estimate the asymmetric drift (Binney & Tremaine 2008): since V²_c = R ∂Φ/∂R, we find

$$\begin{equation} V_c^2 - \bar{V}^2 = \sigma _V^2 -\sigma _U^2\,\left[1 + \frac{\partial \ln \left(\nu \sigma _U^2\right)}{\partial \ln R} + \frac{1}{\nu }\,\frac{R}{\sigma _U^2}\,\frac{\partial \left(\nu \sigma ^2_{UW}\right)}{\partial Z}\right]. \vspace*{6pt} \end{equation} \tag{ 3 }$$

Now assume that (1) the dispersions squared (σ²_U, σ_V², σ²_UW) decline exponentially with radius with scale length h_σ and (2) the tracer density is an exponential function of radius and height with scale lengths h_R and h_Z, that is, ν(R, Z)∝exp (− R/h_R − |Z|/h_Z). The second assumption is only accurate at heights above a few hundred pc; closer to the midplane the exponential form is not accurate, but this is not a concern since we are interested in the region |Z| ≳ 1 kpc. These assumptions were also made by MB12, and in addition they assumed that h_σ = h_R. Equation (3) then becomes

$$\begin{eqnarray} V_c^2 - \bar{V}^2 &=& \sigma _V^2\,{+}\,\sigma _U^2\,\left[ R\left(\frac{1}{h_R}+\frac{1}{h_\sigma }\right) -1\right]\nonumber\\ && + \frac{R}{h_Z}\,\sigma ^2_{UW} -R\,\frac{\partial \sigma ^2_{UW}}{\partial Z}. \end{eqnarray} \tag{ 4 }$$

To evaluate this, we use the expressions for the moments of the velocity distribution from Moni Bidin et al. (2012a),

$$\begin{equation} \sigma _U(R_0,Z) = (82.9 \pm 3.2) + (6.3 \pm 1.1)\cdot (|Z|/\mathrm{kpc}-2.5)\ \mathrm{km\ s}^{-1} \end{equation} \tag{ 5 }$$

$$\begin{equation} \sigma _V(R_0,Z) = (62.2 \pm 3.1) + (4.1 \pm 1.0)\cdot (|Z|/\mathrm{kpc}-2.5)\ \mathrm{km\ s}^{-1} \end{equation} \tag{ 6 }$$

$$\begin{equation} \sigma _W(R_0,Z) = (40.6 \pm 0.8) + (2.7 \pm 0.3)\cdot (|Z|/\mathrm{kpc}-2.5)\ \mathrm{km\ s}^{-1}, \end{equation} \tag{ 7 }$$

and that for σ²_UW from MB12,

$$\begin{eqnarray} \sigma ^2_{UW}(R_0,Z) &=& (1522 \pm 100)\nonumber\\ && + (366 \pm 30)\cdot (|Z|/\mathrm{kpc} - 2.5)\ \mathrm{km}^2\ \mathrm{s}^{-2}. \quad\quad \end{eqnarray} \tag{ 8 }$$

We also take h_z = 0.9 kpc and h_R = h_σ = 3.8 kpc as in MB12, although our own work suggests somewhat smaller values h_R = 2 kpc and h_σ = 3.5 kpc (Bovy et al. 2012b). We find that the solution to Equation (4) can be fit by the formula

$$\begin{equation} V_c^2-\bar{V}^2 = (191\,\mbox{km s}^{-1})^2\,[1+0.19\,(|Z|/\mathrm{kpc}-2.5)], \end{equation} \tag{ 9 }$$

with an rms error of less than 2% for |Z| = 1–4 kpc. Equation (9) implies $V_c-\bar{V}(1.5\ \mathrm{kpc}) = 80$ km s⁻¹ and $V_c-\bar{V}(3.5\ \mathrm{kpc}) = 150\,\mbox{km s}^{-1}$, in good agreement with the corresponding measured values in Moni Bidin et al. (2012a), 70 ± 13 km s⁻¹ and 130 ± 16 km s⁻¹.

We can now ask what Equation (4) predicts for the radial variation of $\bar{V}$. Taking the radial derivative of Equation (4), we find that

$$\begin{eqnarray} 2\,V_c\,\frac{\partial V_c}{\partial R} - 2\,\bar{V}\,\frac{\partial \bar{V}}{\partial R} &=& {-}\frac{\big (V_c^2 - \bar{V}^2 \big)}{h_\sigma } + \sigma _U^2\,\left(\frac{1}{h_R} + \frac{1}{h_\sigma }\right)\nonumber\\ && + \frac{1}{h_Z}\,\sigma ^2_{UW}-\frac{\partial \sigma ^2_{UW}}{\partial Z}\,\nonumber \\ &=& \big(V_c^2 - \bar{V}^2 \big)\left(\frac{1}{R}-\frac{1}{h_\sigma }\right)+\frac{\sigma _U^2}{R}-\frac{\sigma _V^2}{R}. \nonumber\\ \end{eqnarray} \tag{ 10 }$$

We can then estimate $\partial \bar{V}/\partial R$ for a flat circular-speed curve, by setting ∂V_c/∂R = 0:

$$\begin{eqnarray} \bar{V}\,\frac{\partial \bar{V}}{\partial R} &=& 110\,\mbox{km s}^{-1}\times 21\,\mbox{km s}^{-1}\,\mbox{kpc}^{-1}\, \nonumber\\ && \times [1+0.2\,(|Z|/\mathrm{kpc}-2.5)], \end{eqnarray} \tag{ 11 }$$

with an rms error of less than 1% for |Z| = 1–4 kpc. In this equation, 110 km s⁻¹ and 21 km s⁻¹ kpc⁻¹ are $\bar{V}$ and $\partial \bar{V}/\partial R$ at |Z| = 2.5 kpc. We find that $\partial \bar{V}/\partial R$ is 7 km s⁻¹ kpc⁻¹ at Z = 0, growing to 11 km s⁻¹ kpc⁻¹ at |Z| = 1 kpc and 40 km s⁻¹ kpc⁻¹ at |Z| = 3.5 kpc. If we use h_σ = 3.5 kpc (Bovy et al. 2012b), the gradients are larger by about 20%.

For comparison, MB12 estimate that a gradient $\partial \bar{V}/\partial R=17\,\mbox{km s}^{-1}\,\mbox{kpc}^{-1}$ is needed to make their analysis consistent with the expected amount of dark matter, which is close to our estimate of 21 km s⁻¹ kpc⁻¹. MB12 dismiss this possibility, apparently because they confuse constraints on ∂V_c/∂R with constraints on $\partial \bar{V} / \partial R$, as there are, in fact, no direct observational constraints on $\partial \bar{V} / \partial R$.

While $\partial \bar{V} / \partial R$ is not measured for the MB12 sample, constraints on this gradient do exist for a slightly more metal-rich sample of stars at moderate distances from the midplane, which we can use to show that our estimate of $\partial \bar{V}/\partial R$ using Equation (10) agrees with observations. Casetti-Dinescu et al. (2011) study a sample of metal-rich red clump stars at 1.0 kpc ⩽ |Z| ⩽ 2.5 kpc, selected using 0.6 ⩽ J − K_s ⩽ 0.7, while the MB12 data were selected using 0.7 ⩽ J − K_s ⩽ 1.1 to obtain a lower-metallicity, larger-distance sample. Casetti-Dinescu et al. (2011) report that $(\sigma _U,\sigma _V,\bar{V}) = (60,42,180)\,\mbox{km s}^{-1}$ at |Z| ≈ 1 kpc, such that Equation (10) assuming a flat circular-speed curve predicts that $\partial \bar{V}/\partial R = 6.5\,\mbox{km s}^{-1}\,\mbox{kpc}^{-1}$; similarly, they find $(\sigma _U,\sigma _V,\bar{V}) = (80,60,150)\,\mbox{km s}^{-1}$ at |Z| ≈ 1.5 kpc and $(\sigma _U,\sigma _V,\bar{V}) = (90,65,140)\,\mbox{km s}^{-1}$ at |Z| ≈ 2 kpc, such that we predict that $\partial \bar{V}/\partial R = 13\,\mbox{km s}^{-1}\,\mbox{kpc}^{-1}$ and 15 km s⁻¹ kpc, respectively. These predictions are in good agreement with the measurements of Casetti-Dinescu et al. (2011), who find 6.0 ± 1.5 km s⁻¹ kpc⁻¹, 13.0 ± 4.5 km s⁻¹ kpc⁻¹, and 12.0 ± 8.1 km s⁻¹ kpc⁻¹ at these |Z|.

The discussion so far in this section has assumed that V_c is constant with Z and R, whereas observations only show that V_c is independent of R at Z = 0. The circular velocity declines with |Z| for any reasonable local mass distribution, as the following argument shows. As V²_c = −RF_R, we can write

$$\begin{equation} \frac{\partial V_c^2}{\partial Z} = -\frac{\partial }{\partial Z}\,\left({\it RF_R}\right) = -R\frac{\partial F_Z}{\partial R}, \end{equation} \tag{ 12 }$$

where we have used the fact that F_R and F_Z are both derivatives of the potential Φ, such that

$$\begin{equation} \frac{\partial F_R}{\partial Z} = -\frac{\partial }{\partial Z} \left(\frac{\partial \Phi }{\partial R}\right) =-\frac{\partial }{\partial R} \left(\frac{\partial \Phi }{\partial Z}\right) = \frac{\partial F_Z}{\partial R}. \end{equation} \tag{ 13 }$$

As shown below in Section 3, to a good approximation F_Z = −2πGΣ(R, Z) (for Z > 0), where Σ(R, Z) is the surface density. Therefore,

$$\begin{equation} 2\,V_c\,\frac{\partial V_c}{\partial Z} \simeq 2\pi GR\,\frac{\partial \Sigma (R,Z)}{\partial R}=-2\pi G \frac{R}{h_\Sigma }\,\Sigma (R,Z), \end{equation} \tag{ 14 }$$

assuming that Σ(R, Z) declines exponentially with radius with scale length $h_\Sigma$, which we assume is equal to 3.5 kpc (Bovy et al. 2012b). This means that near the midplane, where the density is approximately 0.1 M_☉ pc⁻³,

$$\begin{equation} \frac{\partial V_c}{\partial Z} \simeq -2.8\, \mathrm{km\ s}^{-1}\ \mathrm{kpc}^{-1}\, \left(\frac{Z}{100\,\mathrm{pc}}\right)\quad |Z| \lesssim h_Z, \end{equation} \tag{ 15 }$$

while for |Z| > 1.5 kpc, where Σ(R₀, Z) ≃ 50 M_☉ pc⁻²,

$$\begin{equation} \frac{\partial V_c}{\partial Z} \simeq -7\, \mathrm{km\ s}^{-1}\ \mathrm{kpc}^{-1}\quad |Z| \gg h_Z. \end{equation} \tag{ 16 }$$

Taken together, these results imply that V_c likely does not decrease by more than 30 km s⁻¹ out to |Z| = 4 kpc. Such a decrease does not change the conclusions of this section.

A different way of illustrating the inconsistency of MB12's assumptions is to ask what the radial behavior of V_c has to be such that $\partial \bar{V} / \partial R = 0$. The necessary ∂V_c/∂R decreases from approximately −5 km s⁻¹ kpc⁻¹ at Z = 0, which is consistent with local measurements of the slope of the circular-speed curve, to −14 km s⁻¹ kpc⁻¹ at |Z| = 4 kpc; for h_σ = 3.5 kpc, the gradients are larger in absolute value by 20%. Such a steep drop of the circular velocity with R is about the value of a Keplerian drop-off, −(1/2)V_c/R₀ = −14 km s⁻¹ kpc⁻¹. Such steep gradients are inconsistent with observational evidence: (1) in the limiting case where there is no dark matter halo and all the mass is in an exponential disk (with parameters given in Figure 1), |∂V_c/∂R| is only 3 km s⁻¹ kpc⁻¹ at |Z| = 4 kpc, and if a dark halo is present, the value is even smaller; (2) we show below in Section 3 that ∂V_c/∂Rincreases as |Z| grows, whereas maintaining constant $\bar{V}$ requires that ∂V_c/∂R decreases with |Z|.

Zoom In Zoom Out Reset image size

Having shown that the assumption of a radially constant $\bar{V}$ is suspect, we now ask what the MB12 data do have to say about the surface density at large height |Z|. Following MB12, we start from the Poisson equation in cylindrical coordinates,

$$\begin{equation} \Sigma (R,Z) = -\frac{1}{2\pi G}\,\left[\int _{0}^Z dz\,\frac{1}{R} \,\frac{\partial (RF_R)}{\partial R} + F_Z(R,Z)\right], \end{equation} \tag{ 17 }$$

where we have assumed symmetry around the Z = 0 plane. F_Z(R, Z) is obtained from the steady-state vertical Jeans equation

$$\begin{equation} F_Z(R,Z) = -\frac{\partial \Phi (R,Z)}{\partial Z} = \frac{1}{\nu }\,\frac{\partial \big(\nu \sigma _W^2 \big)}{\partial Z} + \frac{1}{R\nu }\,\frac{\partial \left(R\nu \sigma ^2_{UW}\right)}{\partial R}. \end{equation} \tag{ 18 }$$

One can then proceed (e.g., Kuijken & Gilmore 1989) by approximating the integrand in the first term in square brackets in Equation (17) by its value in the plane. This is zero for a flat circular-speed curve, as RF_R = −V²_c. Kuijken & Gilmore (1989) show that this approximation is good to a few percent at |Z| ≲ 1.5 kpc for any plausible mass distribution. We now revisit this question at larger heights.

Figure 1 shows the error introduced by approximating ∂(RF_R)/∂R by its value at Z = 0, for three different mass distributions: an exponential disk with parameters that are representative of the Milky Way's disk (Bovy et al. 2012a, 2012c), a Navarro–Frenk–White (NFW; Navarro et al. 1997) halo, and a combination of the two in which the disk provides 85% of the radial force at R₀ (this yields a flat circular-speed curve at R₀). This figure shows the fractional difference between the true surface density and that calculated from Equation (17) by approximating ∂(RF_R)/∂R by its value in the plane. The errors introduced by this approximation can be as large as 15% at 4 kpc above the plane. Thus, a systematic uncertainty of 10%–20% in Σ(4 kpc) is unavoidable without fully modeling the potential.

It is, however, possible to show that neglecting the Z dependence of ∂V_c/∂R almost always leads to an underestimate of the surface density, so this approximation gives a robust lower limit on the surface density. We start by writing the Z derivative of the integrand in Equation (17) as

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{1}{R}\,\frac{\partial (R F_R)}{\partial R}\right) = \frac{\partial }{\partial Z} \left(\frac{F_R}{R} + \frac{\partial F_R}{\partial R}\right). \end{equation} \tag{ 19 }$$

If we again use the fact that ∂F_R/∂Z = ∂F_Z/∂R, we find that

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{1}{R}\,\frac{\partial (R F_R)}{\partial R}\right) = \frac{1}{R}\,\frac{\partial F_Z}{\partial R} + \frac{\partial ^2 F_Z}{\partial R^2}. \end{equation} \tag{ 20 }$$

If the first term in square brackets in Equation (17) is neglected, we have F_Z = −2πGΣ(R, Z) (at Z > 0). This approximation is quite accurate if there is substantial mass in a thin disk and |Z| is small, or if the circular-speed curve is nearly flat, since then RF_R is almost independent of radius (see Figure 1). If we then drop this term and approximate the radial derivatives as those of the exponential disk, Equation (20) becomes

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{1}{R}\,\frac{\partial (R F_R)}{\partial R}\right) \simeq \frac{2\pi G}{R\,h_\Sigma } \,\left(1 - \frac{R}{h_\Sigma }\right)\,\Sigma (R,Z). \end{equation} \tag{ 21 }$$

Since $R_0/h_\Sigma$ is 2 or more in any reasonable model for the Milky Way, the right side of this equation is negative, which in turn implies that the integrand in Equation (17) decreases with increasing |Z|. Thus, replacing F_R by its value in the midplane leads to an underestimate of the surface density, and dropping the integral in Equation (17) when the circular-speed curve is locally flat also leads to an underestimate of the surface density.

We can further work out Equation (21) to find the change in the radial slope of the circular-speed curve with height

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{\partial V_c}{\partial R}\right) \simeq \frac{\pi G}{h_\Sigma \,V_c}\,\left(\frac{R}{h_\Sigma }-1\right) \Sigma (R,Z), \end{equation} \tag{ 22 }$$

a result that can also be derived by differentiating Equation (14). Using the same reasoning that led to Equations (15) and (16), we find

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{\partial V_c}{\partial R}\right) \simeq 0.45\,\mathrm{km\ s}^{-1}\ \mathrm{kpc}^{-2}\, \left(\frac{Z}{100\,\mathrm{pc}}\right),\,\quad |Z| \lesssim h_Z, \end{equation} \tag{ 23 }$$

and

$$\begin{equation} \frac{\partial }{\partial Z} \left(\frac{\partial V_c}{\partial R}\right) \simeq 1.1\,\mathrm{km\ s}^{-1}\ \mathrm{kpc}^{-2},\quad |Z| \gg h_Z. \end{equation} \tag{ 24 }$$

Thus, the radial gradient of the circular speed remains close to its value in the plane throughout the region |Z| < 4 kpc that we are investigating.

We have shown that approximating the integrand in Equation (17) as constant with Z leads to a robust lower limit on the surface density as a function of Z. Armed with this result, we can revisit the MB12 analysis. Since the circular-speed curve is flat in the midplane, the integrand is zero in the midplane, so the approximation that the integrand is constant implies that the integral can be neglected. Substituting Equation (18) into Equation (17), we find

$$\begin{eqnarray} &&\Sigma (Z) = -\frac{1}{2\pi G}\,\left[\! -\frac{1}{h_Z}\,\sigma _W^2 + \frac{\partial \sigma _W^2}{\partial Z} + \sigma ^2_{UW}\,\left(\! \frac{1}{R} - \frac{1}{h_R} - \frac{1}{h_\sigma } \!\right) \!\right], \vspace*{6pt}\nonumber\\ \end{eqnarray} \tag{ 25 }$$

where we distinguish between the radial scale length h_R of the tracer population and the radial scale length h_σ of the dispersion squared. MB12 assumed h_R = h_σ = 3.8 kpc, but recent direct measurements of the thicker populations of stars in the Milky Way have shown that their radial scale length is much shorter than that of the dominant, thinner components, such that h_R = 2 kpc is more accurate; we also prefer a shorter scale length for the dispersions, h_σ = 3.5 kpc (Bovy et al. 2012b, 2012c). We will nevertheless show results mostly for the values h_R = h_σ = 3.8 kpc assumed by MB12. We also follow MB12 in assuming h_Z = 900 pc (see below). Note that the surface density as calculated using Equation (25) does not depend on the value of the local circular speed V_c, and our estimate of the local dark matter density is therefore not affected by the uncertainty in V_c.

Using the kinematics from Moni Bidin et al. (2012a) and MB12 in Equation (25), we get the curve labeled as "correct ∂V_c/∂R = 0 approximation" in Figure 2. As we have discussed, this gives a lower limit to the surface density. The gray shaded region in the figure shows the range of surface densities obtained after including the radial integral in Equation (17), using the bottom and top curves in Figure 1. These can be compared to the curve labeled "incorrect $\partial \bar{V} / \partial R = 0$ approximation," which is the curve presented in MB12 as their primary result (their Figure 1).

Zoom In Zoom Out Reset image size

We can also compare these curves with Galactic mass models in the literature, showing the same mass models as MB12. The estimated contribution to Σ(Z) from baryonic matter (VIS for "visual" in Figure 2) is composed of a thin 13 M_☉ pc⁻² layer of interstellar medium, plus a stellar halo and "thick" and "thin" disk components with parameters taken from Jurić et al. (2008); in this model the normalization is chosen so that Σ_disk(1.1 kpc) = 40 M_☉ pc⁻² (Holmberg & Flynn 2004). Figure 2 also shows the expected surface density when various dark matter halo models are added to this baryonic mass. The OM halo model has a profile

$$\begin{equation} \rho _{\mathrm{DM}}(R_0,Z) = \rho _c\left(\frac{R_c^2}{R^2_c + R_0^2 + Z^2}\right), \end{equation} \tag{ 26 }$$

with R_c = 8.01 kpc and ρ_c = 0.0103 M_☉ pc⁻³ (Olling & Merrifield 2001), such that ρ_{☉, DM} = 0.0084 M_☉ pc⁻³. The other halo models were all taken from Weber & de Boer (2010) and have densities following

$$\begin{equation} \rho _{\mathrm{DM}}(r) = \rho _{\odot,{\rm DM}}\,\left(\frac{r}{R_0}\right)^{-\alpha }\,\left(\frac{1 + (r/R_c)^\beta }{1+(R_0/R_c)^\beta }\right)^{-\gamma }, \end{equation} \tag{ 27 }$$

where r is the Galactocentric spherical radius. The standard halo model (SHM) has an NFW profile (α = β = 1, γ = 2) with R_c = 10.8 kpc and ρ_{☉, DM} = 0.0084 M_☉ pc⁻³. The N97 model also has an NFW profile with R_c = 20 kpc and ρ_{☉, DM} = 0.0061 M_☉ pc⁻³. The MIN model has a profile as in Equation (26) (i.e., α = 0, β = 2, γ = 1) with R_c = 5 kpc and ρ_c = 0.019 M_☉ pc⁻³ (Weber & de Boer 2010), such that ρ_{☉, DM} = 0.0053 M_☉ pc⁻³. All of these models were constrained by assuming that V_c(R₀) = 244 ± 10 km s⁻¹, except for the OM model, which has V_c(R₀) = 220 km s⁻¹. The difference between the predictions of these models for Σ(Z) in Figure 2 is mainly due to their different value for ρ_{☉, DM}.

It is clear from Figure 2 that the MB12 data are fully consistent with the predictions from several standard dark matter models when using the correct assumption ∂V_c/∂R = 0. In particular, the values of the surface density throughout the range |Z| = 1–4 kpc are consistent with the SHM, which has ρ_DM ≃ 0.01 M_☉ pc⁻³ = 0.38 GeV cm⁻³. Figure 2 also shows the effect of using the more appropriate radial scale length h_R = 2 kpc, which increases the surface density even further. The slope of the measured Σ(Z) implies a minimum dark matter density of ρ_DM = 0.007 ± 0.001 M_☉ pc⁻³ = 0.27 ± 0.04 GeV cm⁻³ for h_R = 3.8 kpc and ρ_DM = 0.0085 ± 0.0015 M_☉ pc⁻³ = 0.32 ± 0.06 GeV cm⁻³ for h_R = 2.0 kpc. These uncertainties are statistical and do not include the systematic uncertainty associated with the approximation that the circular-velocity curve is flat at all heights, which adds at most 0.1 GeV cm⁻³ (see Figure 1), or that associated with the measurement of the dispersions (see, e.g., the erratic behavior of σ²_UW in Figure 8 of Moni Bidin et al. 2012a), or the assumption that the tracers follow an exponential distribution in radius and height with the assumed scale lengths. These systematic uncertainties are at least of the same magnitude as the statistical uncertainties. For example, if instead of h_Z = 900 pc we use h_Z = 700 pc, which is the mass-weighted mean scale height in the range 2–4 kpc using the measurements of Bovy et al. (2012a, 2012c), ρ_DM increases by 0.001 M_☉ pc⁻³; using a varying mass-weighted h_Z as a function of |Z| as given by those same measurements changes ρ_DM by 0.002 M_☉ pc⁻³. Taking these systematics into account, our best estimate for the local dark-halo density is ρ_DM = 0.008 ± 0.003 M_☉ pc⁻³ = 0.3 ± 0.1 GeV  cm⁻³.

In this paper, we have shown that the assumption of a radially constant mean azimuthal velocity for the stellar tracers used by MB12 is physically implausible. This assumption is the reason why the MB12 analysis finds a constant Σ(Z) at 2 kpc ≲ |Z| ≲ 4 kpc, in seeming contradiction with the standard expectation of a dark matter density ρ_DM ≈ 0.01 M_☉ pc⁻³ = 0.38 GeV cm⁻³ in this region. The mean azimuthal velocity $\bar{V}$ of the warm tracer population of stars used by MB12—stars that are part of the thicker components of the Milky Way disk that reach large heights above and below the plane—is significantly different from the circular velocity V_c, with $\bar{V} \lesssim 0.5V_c$ at |Z| = 3 kpc. We have shown in Section 2 that for a circular-speed curve that is close to flat, the expected $\partial \bar{V} / \partial R$ reaches tens of  km s⁻¹ kpc⁻¹ at a few kpc from the Galactic midplane. Indeed, the assumption that $\partial \bar{V} / \partial R = 0$ requires that the circular-speed curve falls off in a Keplerian manner at a few kpc above the plane, in clear contradiction with observations.

We derived an alternative formula for Σ(Z) that assumes that the circular-velocity curve is flat at all heights above the plane, and we showed that this approximation leads to a lower limit for all plausible mass distributions. This approximation sidesteps the issue of the unknown radial trend of $\bar{V}$ and as such makes no assumptions about it. Applying this alternative formula, we find that the MB12 data give a lower limit that is fully consistent with the standard local density of dark matter of ρ_DM ≈ 0.01 M_☉ pc⁻³, and that they imply a local dark matter density of 0.008 ± 0.003 M_☉ pc⁻³ = 0.3 ± 0.1 GeV  cm⁻³, where the error bars include both statistical and less well known systematic errors. Therefore, our analysis shows that the locally measured density of dark matter is consistent with that extrapolated from halo models constrained at Galactocentric distances r ≳ 20 kpc: for example, extrapolating the best-fit halo profiles of Xue et al. (2008) and Deason et al. (2012), obtained from fitting halo stars with 20 kpc ≲ r ≲ 60 kpc, gives ρ_DM ≈ 0.006–0.011 M_☉ pc⁻³.

The breakdown of the assumptions made in this simple, "model-independent" Jeans analysis is such that the measurement has a systematic uncertainty reaching 10%–20% at |Z| = 4 kpc. Therefore, a precise determination of the local dark matter density from observations at large Z using the Jeans analysis of MB12 requires data that span a wide range in R such that the radial gradient of the velocity moments, foremost $\bar{V}$, can be determined. The Gaia mission (Perryman et al. 2001) will provide such measurements in the near future. Currently, data that span multiple kpc in R and Z are available from the SDSS/SEGUE survey (Yanny et al. 2009). In contrast to the data from Moni Bidin et al. (2012a), the SEGUE volume selection is known in detail (e.g., Bovy et al. 2012c), so both the spatial structure and kinematics can be obtained for the same set of tracer stars. The SEGUE data also have the advantage that they contain information on the elemental abundances of each star, allowing the tracers to be divided into sets of stars with simple distribution functions (e.g., Bovy et al. 2012b, 2012c).

Finally, we note that our estimate is of the mean halo density between 1 and 4 kpc above the Galactic midplane. The halo density at the Sun, which is the relevant quantity for direct dark matter detection experiments, is likely to be larger because of two effects. The density in the midplane for a spherical NFW halo with a scale radius of 22.25 kpc (Xue et al. 2008) is 7% larger than at |Z| = 2.5 kpc. Besides this purely geometric effect, the gravitational influence of the disk further increases the midplane dark matter density. An isothermal halo with isotropic velocity dispersion σ has a density $\rho \propto \int d\vec{v}\, e^{-E/\sigma ^2} = e^{-\Phi /\sigma ^2}$, so we expect that

$$\begin{eqnarray} \frac{\rho _{{\rm DM}}(|Z|)}{\rho _{{\rm DM}}(0)} - 1 &\simeq & {-}\frac{2\pi G\Sigma _{\rm disk}|Z|}{\sigma ^2}\nonumber\\ &=& \,{-}0.20\,\frac{\Sigma _{\rm disk}}{50\,M_\odot \,\mbox{pc}^{-2}} \left(\! \frac{130\,\mbox{km s}^{-1}}{\sigma } \!\right)^2\frac{|Z|}{2.5\,\mbox{kpc}}. \nonumber\\ \end{eqnarray} \tag{ 28 }$$

These two effects imply that the dark matter density in the midplane is enhanced over the value derived in this paper by about 30%. This agrees with the estimated enhancement in an N-body simulation by Garbari et al. (2011).

It is a pleasure to thank the anonymous referee, Dana Casetti-Dinescu, Christian Moni Bidin, Jerry Ostriker, Justin Read, Hans-Walter Rix, Martin Smith, and Nadia Zakamska for helpful comments. Support for Program number HST-HF-51285.01-A was provided by NASA through a Hubble Fellowship grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555.