KINEMATICAL AND CHEMICAL VERTICAL STRUCTURE OF THE GALACTIC THICK DISK. II. A LACK OF DARK MATTER IN THE SOLAR NEIGHBORHOOD^*,†

Today, it is widely accepted that a dominant fraction of the mass in the universe is in the form of a non-luminous (dark) matter (DM), whose nature is still unknown. Despite this general agreement, and decades of investigations, its spatial distribution and its density in the solar neighborhood (ρ_{☉, DM}) are still poorly constrained by the observations. The local density⁵ extrapolated from the Milky Way rotation curve and other observational data, assuming a spherical DM halo, is in the range 5–13 mM_☉ pc⁻³ (Weber & de Boer 2010). However, ρ_{☉, DM} could be significantly different in the case of an oblate or prolate halo, or in the presence of a DM disk. Our poor knowledge of the DM density distribution is very unfortunate, since this information is crucial to clarify the nature and properties of DM itself. The shape of the dark halo, quantified by its shortest-to-longest axis ratio q, is in fact an important diagnostic on the nature of the DM: for example, a round halo (q ≈ 1) is expected by hot DM models (Peebles 1993), while a very flat halo (q ≈ 0.2) is preferred if cold molecular gas or massive decaying neutrinos are the main constituent of the DM (Pfenniger et al. 1994; Sciama 1990). Noticeably, great effort is currently spent in experiments for direct DM detection. The results of these experiments are degenerate between the unknown interaction cross section of the searched particles and their local density (e.g., Gaitskell 2004; Aprile et al. 2005). Thus, most of the works estimating the properties of the weakly interacting massive particles have so far drawn their conclusions assuming the local DM density of the standard halo model (SHM; Jungman et al. 1996), ρ_{☉, DM} = 8 mM_☉ pc⁻³. Clearly, observational constraints on the DM density at the solar position and on the flattening of the dark halo are key components for revealing the secrets of the Galactic DM.

"Weighing" the Galactic disk by means of the spatial distribution and kinematics of its stellar component is an ancient art, dating back nearly one century (Kapteyn 1922; Oort 1932). The difference between the measured mass and the visible mass provides an estimate of the amount of DM in the volume under analysis, and constraints on the shape of the DM halo can be derived. The fundamental basis for this classical measurement is the application of the Poisson–Boltzmann and Jeans equations to a virialized system in steady state. This allows us to estimate either the local density at the solar position or the surface density (mass per unit area) of the mass within a given volume. Garrido Pestaña & Eckhardt (2010) recently argued against the reliability of this method, but Sánchez-Salcedo et al. (2011) demonstrated that it can be applied to the Galactic disk because it is in equilibrium with the Galactic potential. The measurements of the dynamical mass in the solar neighborhood are abundant in the literature, and all but a few works (Bahcall 1984; Bahcall et al. 1992) came to the same overall conclusion that "there is no evidence for a significant amount of disk DM." Different interpretations of this statement have been presented, however: while the measurements of the surface density usually match the expectations of the visible mass plus a classical spherical DM halo (Kuijken & Gilmore 1989a; Flynn & Fuchs 1994; Siebert et al. 2003; Holmberg & Flynn 2004; Bienaymé et al. 2006), the estimates of the local volume density usually find a much lower quantity of DM (Kuijken & Gilmore 1989b; Crézé et al. 1998; Holmberg & Flynn 2000; Korchagin et al. 2003; de Jong et al. 2010), although they are compatible with the presence of a classical halo within the uncertainties. Despite many decades of investigation, little progress has been made beyond this statement: for example, the most recent measurements are still compatible with both the expectations of the SHM and the complete absence of DM in the solar neighborhood (e.g., Garbari et al. 2011). Some weak constraints on the properties of the DM halo were given by Crézé et al. (1998) and Bienaymé et al. (2006), who claimed that a DM halo flattening q ⩽ 0.5 is excluded by the observations.

The main limitations to the measurements of the dynamical mass are imposed by the great observational effort required, because the information about the kinematics and spatial distribution of a vast number of stars is needed. Some approximations were thus always introduced in the formulation, accurate only at low Galactic heights, whose validity has since been questioned in the literature. For example, Siebert et al. (2008) and Smith et al. (2009) found that at only Z = 1 kpc from the Galactic plane the potential is not separable in the radial and vertical coordinates, and neglecting the non-diagonal term of the dispersion matrix could introduce a bias. Garbari et al. (2011) found that this becomes an issue at only Z = 0.5 kpc and argued against the use of the approximate formulation of Holmberg & Flynn (2000, 2004) at high Z. As a consequence, all previous investigations were limited to Z ⩽ 1.1 kpc from the Galactic plane. The amount of DM expected in this volume is small compared to observational uncertainties, and no strong conclusion could be drawn. Moreover, only Korchagin et al. (2003) directly calculated the mass density from an analytical expression, while the other investigations estimated it by comparison of the observational quantities with the expectations of a Galactic mass model.

In this paper, we propose a new approach to estimate the surface mass density of the Galactic disk. Our formulation leads to an analytical expression that is exact within the limits of validity of a few basic assumptions. The surface density can thus be directly calculated at any distance from the Galactic plane. The calculation requires a knowledge of the spatial variations of the three-dimensional kinematics of a test stellar population. Our calculations are based on the results of our recent investigation of the Galactic thick disk kinematics (Carraro et al. 2005; Moni Bidin 2009). Similar results were presented by Moni Bidin et al. (2010), who searched for a signature of the DM disk, predicted by the merging scenario of thick disk formation.

In the following, we will use the cylindrical Galactic coordinates (R, θ, Z), where R is the Galactocentric distance, θ is directed in the sense of Galactic rotation, and Z is positive toward the north Galactic pole. The respective velocity components are ($\dot{R},\dot{\theta },\dot{Z}$) = (U, V, W).

Our formulation is based on the integrated Poisson equation in cylindrical coordinates:

$$\begin{equation} {-}4\pi G\Sigma (Z) = \int _{-Z}^{Z} \frac{1}{R}\frac{\partial }{\partial R}(R F_{R}) dz +2\cdot [F_{z}(Z) - F_{z}(0)], \end{equation} \tag{ 1 }$$

where G is the gravitational constant, Σ(Z) is the surface density of the mass comprised between ±Z from the Galactic plane, and the radial and vertical component of the force per unit mass, F_R and F_Z, respectively, can be expressed through the Jeans equations

$$\begin{eqnarray} F_{R}&=& {-}\frac{\partial \phi }{\partial R}=\frac{1}{\rho }\frac{\partial (\rho \overline{U^2})}{\partial R}+ \frac{1}{\rho }\frac{\partial (\rho \overline{UW})}{\partial Z}+\frac{\overline{U^2}-\overline{V^2}}{R} \nonumber\\ &&+ \frac{1}{\rho }\frac{\partial (\rho \overline{U})}{\partial t}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&F_{Z}=-\frac{\partial \phi }{\partial Z}=\frac{1}{\rho } \left[\frac{\partial (\rho \overline{W^2})}{\partial Z}+ \frac{\rho \overline{UW}}{R}+\frac{\partial (\rho \overline{UW})}{\partial R}+ \frac{\partial (\rho \overline{W})}{\partial t} \right], \nonumber\\ \end{eqnarray} \tag{ 3 }$$

where ρ is the volume density and ϕ is the gravitational potential. The following basic assumptions can be made, for symmetry reasons, for the stellar population under study.

• I.  
  The test population is a virialized system in steady state. All the temporal derivatives are therefore null:

  $$\begin{equation} \frac{\partial (\rho \overline{U})}{\partial t}=\frac{\partial (\rho \overline{W})}{\partial t}=0. \end{equation} \tag{ 4 }$$
• II.  
  The vertical component of the force per unit mass is null on the plane:

  $$\begin{equation} F_{Z}(0)=0. \end{equation} \tag{ 5 }$$
• III.  
  The trend of velocity dispersions with Z is symmetric with respect to the plane:

  $$\begin{equation} \sigma _i(Z)=\sigma _i(-Z), \end{equation} \tag{ 6 }$$

  with i = U, V, W.
• IV.  
  $\overline{UW}$ is antisymmetric with respect to the Galactic plane:

  $$\begin{equation} \overline{UW}(Z)=-\overline{UW}(-Z). \end{equation} \tag{ 7 }$$
• V.  
  The net radial and vertical bulk motion of the trace population is null:

  $$\begin{equation} \overline{U}=\overline{W}=0. \end{equation} \tag{ 8 }$$

Sánchez-Salcedo et al. (2011) demonstrated the validity of assumption (I) that was recently questioned in the literature (Garrido Pestaña & Eckhardt 2010). The hypotheses (II)–(IV) are required for symmetry reasons, and they imply that the integral of the velocity dispersions in ±Z is twice the integral between zero and Z, while the integral of $\overline{UW}$ is null, and $\overline{UW}$(0) = 0. It can be easily seen that assumption (IV) directly results from the Jeans equations, if the symmetry requirements F_R(Z) = F_R(− Z) and F_Z(Z) = −F_Z(− Z) are to be satisfied. Moni Bidin et al. (2010) also showed that the cross term must be assumed antisymmetric with respect to Z, or else the calculation leads to unphysical results. Recent investigations detected a non-null mean radial motion of stars (Casetti-Dinescu et al. 2011; Moni Bidin et al. 2012). However, this is not a consequence of a global motion of disk stars, but it is rather a signature of local kinematical substructures, so that assumption (V) is still valid as a general property of disk stars.

Inserting Equations (2) and (3) into Equation (1), and making use of the assumptions (I)–(V), the resulting equation can be analytically solved for Σ(Z) at any R and Z, if the three-dimensional analytical expressions for the kinematics and the mass density distribution of the test population are given. To provide this extensive input, we make use of observational results only. The spatial distribution of the test population is fixed by the following assumptions.

• VI.  
  The volume density of the test population decays exponentially both in the radial and vertical direction, with exponential scale height and length h_{Z, ρ} and h_{R, ρ}, respectively:

  $$\begin{equation} \rho (R,z)=\rho _\odot \cdot \exp \left(-\frac{Z-Z_\odot }{h_{Z,\rho }}-\frac{R-R_\odot }{h_{R,\rho }}\right), \end{equation} \tag{ 9 }$$

  where ρ_☉ is the density at (R, Z) = (R_☉, Z_☉).
• VII.  
  The scale length of the test population, h_{R, ρ}, is invariant with respect to Galactic height:

  $$\begin{equation} \frac{\partial h_{R,\rho }}{\partial Z}=0. \end{equation} \tag{ 10 }$$

These hypotheses are a classical representation of Galactic disk-like populations, and (VII) is confirmed by the empirical results of Cabrera-Lavers et al. (2007). Theoretical considerations suggest that the vertical density profile should be much closer to a sech² function than an exponential (Camm 1950, 1952), but observational evidences confirm that the Galactic thick disk, the object of our investigation, is well described by the double exponential law of Equation (9) (e.g., Hammersley et al. 1999; Siegel et al. 2002; Jurić et al. 2008). The assumed vertical decay is surely an accurate fit of the disk density at our large Galactic heights, where even the sech² profile approximates to an exponential decline.

The calculation is simplified by two additional working hypotheses.

• VIII.  
  The rotation curve is locally flat in the volume under study:

  $$\begin{equation} \frac{\partial \overline{V}}{\partial R}=0. \end{equation} \tag{ 11 }$$
• IX.  
  The disk test population is not flared in the volume of interest:

  $$\begin{equation} \frac{\partial h_{Z,\rho }}{\partial R}=0. \end{equation} \tag{ 12 }$$

These are very common assumptions, whose function in our approach is neglecting few small terms, which complicate the formulation while introducing only second-order corrections. The observational evidences in their support will be discussed in Sections 4.2 and 4.4, where the influence of their breakdown will also be analyzed.

It is easily demonstrated that, by means of our assumptions, the final expression for the surface density reduces to

$$\begin{eqnarray} {-}2\pi G\Sigma (R,Z) &=& \frac{\partial \sigma ^2_W}{\partial Z}-\frac{\sigma ^2_W}{h_{Z,\rho }} - \int _{0}^{Z}\frac{\sigma ^2_U}{Rh_{R,\rho }} dz \nonumber\\ &&+\! \int _{0}^{Z}\!\frac{\partial ^2 \sigma ^2_U}{\partial R^2} dz +\! \int _{0}^{Z}\left(\frac{2}{R}\!-\!\frac{1}{h_{R,\rho }}\right)\!\cdot \frac{\partial \sigma ^2_U}{\partial R} dz \nonumber\\ &&- \frac{1}{R}\int _{0}^{Z}\frac{\partial \sigma ^2_V}{\partial R} dz + \overline{UW}\left(\frac{2}{R}-\frac{1}{h_{R,\rho }}\right) \nonumber\\ &&+ 2\frac{\partial \overline{UW}}{\partial R} - \frac{1}{h_z\rho }\int _{0}^{Z}\frac{\partial \overline{UW}}{\partial R} dz. \end{eqnarray} \tag{ 13 }$$

This equation is exact within the limits of validity of the underlying assumptions, and it can be used to calculate the surface density of the mass within ±Z from the Galactic plane at the Galactocentric distance R. The estimate requires the knowledge of the two parameters h_{Z, ρ} and h_{R, ρ}, the vertical trend of the three dispersions and of $\overline{UW}$, and the radial derivative of σ_U(Z), σ_V(Z), and $\overline{UW}(Z)$. A great observational effort, impossible until the last decade, is required to collect this quantity of information. Nevertheless, this will be eventually gathered by modern extensive surveys such as the Sloan Digital Sky Survey (York et al. 2000), Gaia (Wilkinson et al. 2005), and Large Synoptic Survey Telescope (Ivezić et al. 2008). The measurement of the mass distribution in the Galaxy will then be possible with unprecedented detail.

As shown in Section 3, the data used in the present study provide no information about the radial behavior of the kinematics, and we are forced to model it with an additional assumption. We adopt the most likely hypothesis:

• X.  
  $\overline{UW}$ and the square of the dispersions exponentially decay in the radial direction, with the same scale length as that of the volume mass density.

The effects of alternative radial trends will be analyzed in Section 4.1. Assumption (X) is observationally confirmed only for σ_W (van der Kruit & Searle 1981, 1982), and its extension to the other components relies on the assumption of radially constant anisotropy, i.e., the radial constancy of the ratio of the dispersions. However, both theoretical calculations and observations support this hypothesis: Cuddeford & Amendt (1992) demonstrated, by means of numerical integration of orbits, that it is the best representation of the radial trend of dispersions for R  ⩽  9 kpc, and the most detailed observations available to date also confirm it (Lewis & Freeman 1989), although Neese & Yoss (1988) prefer a linear radial decay for σ_U. Moni Bidin et al. (2012) also find that the observational data are consistent with assumption (X).

Using assumption (X), Equation (13) can be simplified substantially, and after a calculation involving only simple integrals and derivatives, we obtain

$${\selectfont\fontsize{9.3}{11.3}{\begin{eqnarray} &&\Sigma (Z)=\nonumber\\ &&\frac{1}{2\pi G}\Bigg[ k_1\cdot \int _{0}^{Z}\sigma _{U}^{2}dz+k_2\cdot \int _{0}^{Z}\sigma _{V}^{2}dz+k_3\cdot \overline{UW} +\frac{\sigma ^{2}_{W}}{h_{Z,\rho }}-\frac{\partial \sigma ^{2}_{W}}{\partial Z}\Bigg],\nonumber\\ \end{eqnarray}}} \tag{ 14 }$$

where

$$\begin{eqnarray} &&k_1=\frac{3}{R_\odot \cdot h_{R,\rho }}-\frac{2}{h_{R,\rho }^{2}}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} &&k_2=-\frac{1}{R_\odot \cdot h_{R,\rho }}, \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} &&k_3=\frac{3}{h_{R,\rho }}-\frac{2}{R_\odot }. \end{eqnarray} \tag{ 17 }$$

Our estimate of the surface mass density is based on the results of Moni Bidin et al. (2012), who measured the kinematics of the Galactic thick disk, and its variation with distance from the plane, between Z = 1.5 and 4.5 kpc. In brief, they studied a sample of ∼400 thick disk red giants toward the south Galactic pole, vertically distributed with respect to the Galactic plane, and derived their three-dimensional kinematics by means of Two Micron All Sky Survey photometry (Skrutskie et al. 2006), SPM3 absolute proper motions (Girard et al. 2004), and radial velocities (Moni Bidin 2009). Thus, the data do not provide any information about the variation of the kinematics with Galactocentric distance, which will be modeled with assumption (X). The surface mass density will therefore be calculated by means of Equation (14). Moni Bidin et al. (2012) detected a clear increment with Z of all the dispersions, well represented by a linear fit. We will adopt for σ_U(Z), σ_V(Z), and σ_W(Z) the relations given by Equations (3)–(5) of Moni Bidin et al. (2012). The vertical trend of $\overline{UW}$ will be taken from the linear fit shown in Figure 2 of Moni Bidin et al. (2010), which yields the solution

$$\begin{equation} \overline{UW}=(1522\pm 100)+(366\pm 30)\cdot (Z-2.5) \,\mathrm{km^2 \,s^{-2}}, \end{equation} \tag{ 18 }$$

where Z is in kiloparsecs.

Measuring the kinematical properties and the spatial distribution of the test population by means of the same observed stars would be highly desirable. Unfortunately, our sample is not suitable for this: even if Girard et al. (2006) used it to estimate the thick disk scale height, their measurement could have been affected by the thin disk contamination unaccounted for in their study (Moni Bidin et al. 2012). However, our sample can be considered representative of the thick disk kinematics at the solar position, whose spatial distribution has been extensively studied in the last three decades. This would not be true if the thick disk was a mixture of sub-populations with different properties. In this case, the spatial distribution of the sub-population under study, or the dominant one, should be preferred. Nevertheless, to our knowledge, there is no evidence in literature supporting such high degree of inhomogeneity in the thick disk, nor that the sample under study is for some reason peculiar. Finally, we fix R = R_☉ = 8.0  ±  0.3 kpc, and h_{R, ρ} = 3.8  ±  0.2 kpc, h_{Z, ρ} = 900  ±  80 pc from the average of 16 and 21 literature measurements, respectively.⁶ The impact of these parameters on the final results will be analyzed in Section 4.1. The errors were defined, as in Moni Bidin et al. (2010), from the statistical error-on-the-mean. If the literature estimates can be considered independent measurements of an underlying quantity, this must be the most rigorous estimate of the true uncertainty. However, the scale height and length of the Galactic thick disk are traditionally poorly constrained, although the measurements converged considerably in the last years. For this reason, Moni Bidin et al. (2010) also considered enhanced errors (0.4 and 0.12 kpc for h_{R, ρ} and h_{Z, ρ}, respectively), possibly more representative of the true uncertainties. In this case, the final errors on Σ(Z) are enhanced by about 50%, and the significance of the results thus decreases by a factor of ∼0.7. Nevertheless, adopting these larger errors would not affect the general conclusions of our work noticeably because, as will appear more clearly later, the significance of the results remains very high even in this case.

The error on Σ(Z) was calculated from the propagation of the uncertainties on the four kinematical quantities and the three parameters. The resulting equation is particularly cumbersome, but it can be obtained through simple derivations of the terms on the right-hand side of Equation (14). Monte Carlo simulations, kindly provided us by the referee, confirmed that the final error is a good estimate of the uncertainty propagated from the ones associated with the quantities entering in the calculations. The simulations were performed repeating the estimate of Σ(Z) after varying the input quantities, each one randomly drawn from a Gaussian distribution with mean and dispersion given by the assumed value and error, respectively.

The resulting vertical profile of the surface mass density Σ(Z) is shown in the upper panel of Figure 1. We find Σ(1.5 kpc) = 55.6 ± 4.7 M_☉ pc⁻², and the profile is nearly flat, increasing by only ∼2 M_☉ pc⁻² up to Z = 3 kpc and eventually bending down, while the error monotonically increases up to 10 M_☉ pc⁻² at Z = 4 kpc. The decrease of Σ(Z) in the last kiloparsec is unphysical, but it is negligible compared to the errors, corresponding to only 0.07σ. This could be corrected by a tiny change of the input parameters, for example, increasing h_{R, ρ} by 0.1 kpc, but we will not manipulate them to avoid the introduction of arbitrariness in the results. If the solution is extrapolated to Z = 0, the unphysical result Σ(0) ≠ 0 is obtained. This is due to the fact that the vertical trend of the kinematical quantities was derived from a linear fit. This is a good approximation in the Z-range under analysis, but the dispersions should depart from it at lower heights, bending down with steeper gradient. The most clear case is $\overline{UW}$, which should be zero on the plane for symmetry reasons, but the linear fit returns $\overline{UW}(0)=607$ km² s⁻². The extrapolation of the results thus fails to account for the changed gradient of the kinematical quantities and, as a consequence, of Σ(Z). In fact, as shown later (Figure 6), a steeper curve of the surface mass density at lower Z (and a lower extrapolation to Z = 0) is recovered when using the steeper gradients of kinematics of Casetti-Dinescu et al. (2011). It can be noted that the linear fit is a worse approximation in this case, as reflected by Σ(Z) being too flat compared to any model, because the dispersions vary their gradient more quickly nearer to the plane. In conclusion, the extrapolation of our results outside the Z-range where they were obtained is not allowed.

Zoom In Zoom Out Reset image size

In Figure 1, the results are compared with the expectations for the known visible mass, indicated by the thick gray curve labeled as VIS. This was modeled as the sum of a thin layer of interstellar medium (ISM) of 13 M_☉ pc⁻² (Holmberg & Flynn 2000), plus three stellar components, the halo, the thin disk, and the thick disk, whose geometrical parameters were taken from Jurić et al. (2008). Their local density at Z = 0 was normalized so that the total surface density of the stellar disk is Σ_disk(1.1 kpc) = 40 M_☉ pc⁻² (Holmberg & Flynn 2004; Bienaymé et al. 2006). This is currently the best estimate assumed in Galactic mass models (e.g., Dehnen & Binney 1998; Olling & Merrifield 2001; Weber & de Boer 2010), but the local mass density of both the ISM and the stellar component is still affected by observational uncertainties. However, as discussed below, the exact value of these parameters does not have a significant influence on our results.

The estimate of the surface mass density matches the expectation of visible mass alone, and the degree of overlap between the two curves is striking. There is no need for any dark component to account for the results: the measured Σ(Z) implies a local DM density ρ_{☉, DM} = 0 ± 1 mM_☉ pc⁻³. This estimate negligibly changes (ρ_{☉, DM} = 0.4 ± 1.2 mM_☉ pc⁻³) if we assume Σ_disk(1.1 kpc) = 35 M_☉ pc⁻² in the visible mass model (Holmberg & Flynn 2000; Garbari et al. 2011). Fitting a DM halo model to the observations is pointless because, disregarding its exact shape, the procedure necessarily converges to a zero-density solution. The more complex task of building a model able to reproduce our results and other observational constraints (e.g., the Galactic rotation curve, the gas disk flare) is beyond the scope of this paper. Nevertheless, we can compare the measurements with the expectations of the most popular DM halo models.

A great quantity of models have been proposed in the literature to describe the spatial distribution of the Galactic DM, but only the mass between 0 and 4 kpc from the plane at R = R_☉ is involved in the comparison with the observations. Hence, the local density ρ_{☉, DM} and the halo flattening q are the only relevant parameters, while the exact functional form is not critical. Indeed, it is easy to see that different models with the same ρ_{☉, DM} and q are indistinguishable in terms of the expected quantity of DM in the volume under analysis. We will therefore compare the observations with a set of selected models, chosen on the basis of the expected ρ_{☉, DM}. We will assume spherical models here (q = 1), while the effects of varying the flattening parameter will be analyzed in more detail in Section 4.7.

Olling & Merrifield (2001) presented a family of self-consistent models, based on a classical non-singular isothermal spheroid:

$$\begin{equation} \rho _\mathrm{DM}(R,Z) = \rho _c\left[\frac{R^2_c}{R^2_c+R^2_\odot +(Z/q)^2}\right] \end{equation} \tag{ 19 }$$

(e.g., van Albada & Sancisi 1986; Kent 1987), where the core radius R_c and the central density ρ_c depend on q, so that the resulting rotation curve is independent of it (Olling 1995). As a conservative choice, we will assume, among the solutions proposed by Olling & Merrifield (2001) with R_☉ = 7.8–8.5 kpc and Σ_disk(1.1 kpc) ⩾ 30 M_☉ pc⁻², the model with the minimum local density (ρ_{☉, DM} ≈ 10 mM_☉ pc⁻³, hereafter model OM), which have R_c = 8.01 kpc and ρ_c = 20.6 mM_☉ pc⁻³.

A more general expression of the DM halo shape is

$$\begin{eqnarray} \rho _\mathrm{DM}(R,Z) &=& \rho _{\odot,{\rm DM}}\cdot \left(\frac{\sqrt{R^2+(Z/q)^2}}{R_\odot }\right)^{-\alpha } \nonumber\\ && \cdot \left[\frac{1+\left(\frac{\sqrt{R^2\,{+}\,(Z/q)^2}}{R_c}\right)^{\beta }}{1+\left(\frac{R_\odot }{R_c}\right)^{\beta }}\right]^{-\gamma }, \end{eqnarray} \tag{ 20 }$$

where the indices (α, β, γ) characterize the radial fall-off of the density distribution. Various sets of indices have been proposed in the literature, suggesting either "cuspy" (e.g., Navarro et al. 1997; Moore et al. 1999; Binney & Evans 2001; Ludlow et al. 2009) or "cored" profiles (e.g., de Boer et al. 2005; Narayan et al. 2005; Gentile et al. 2007; Salucci et al. 2007; Oh et al. 2008). The most representative models were analyzed by Weber & de Boer (2010), who fixed the best-fit parameters on the basis of the most recent observational constraints. Among these, we will consider the model which better fits the observations (in terms of the lowest χ²), the Navarro et al. (1997, NFW) profile (α = β = 1, γ = 2) with R_c = 10.8 kpc and ρ_{☉, DM} = 8.4 mM_☉ pc⁻³. This value coincides with the SHM density usually assumed in direct DM detection experiments (Jungman et al. 1996), and it will be referred to as SHM. We will also consider the NFW profile with the lowest local density (R_c = 20 kpc, ρ_{☉, DM} = 6.1 mM_☉ pc⁻³, hereafter model N97), and the model with the minimum local DM density, a pseudo-isothermal profile (α = 0, β = 2, γ = 1; de Boer et al. 2005, hereafter model MIN) with R_c = 5 kpc and ρ_{☉, DM} = 5.3 mM_☉ pc⁻³. This value is usually assumed as a lower limit for the local DM density (Garbari et al. 2011; Weber & de Boer 2010).

In Figure 1, the measured vertical profile of the surface density is compared to the expectations of the selected DM halo models. The OM model is excluded at the 8σ level, while the SHM at the 6σ level. Even the MIN model, with the minimum local density extrapolated from the Galactic rotation curve (Garbari et al. 2011), is 4.1σ more massive than the observed curve. If we assume Σ_disk(1.1 kpc) = 35 M_☉ pc⁻² in the visible mass model, the results have a lower but still very significant departure: even the MIN model would still be excluded at the 3.6σ level.

It could be argued that integrating the linear vertical profiles of the velocity dispersions in the range Z = 0–4 kpc implicitly assumes their extrapolation down to Z = 0, while they were measured only for Z ⩾ 1.5 kpc. Indeed, the vertical trend of the dispersions is not expected to be strictly linear, but to be shallower at higher distance from the plane (see, for example, the discussion of Casetti-Dinescu et al. 2011). Nevertheless, assuming for Z ⩽ 1.5 kpc the steeper relations measured by Casetti-Dinescu et al. (2011), the resulting surface density decreases by only 0.5 M_☉ pc⁻². Arbitrarily increasing all the gradients by a factor of two for Z = 0–1.5 kpc, Σ(Z) decreases by 5.5 M_☉ pc⁻², enhancing the disagreement with the DM halo models. In conclusion, the integration of the profiles in the entire range is likely a good approximation, and it cannot be the cause of the mismatch between the measurements and the model expectations.

The increment of the surface density between 1.5 and 4 kpc (Σ_>1.5(Z)) can be measured integrating Equations (13) and (14) in this interval, instead of 0–4 kpc. This provides a measurement of the surface density of the mass enclosed between 1.5 and 4 kpc from the plane only. This approach has the advantage of being free of the uncertainties related to the total amount of visible mass and the extrapolation of the kinematics to Z  ⩽  1.5 kpc, although we have already demonstrated that these points do not invalidate our conclusions. In fact, all the ISM and the majority of the stellar mass are found below the range of integration, and a decrease of Σ_disk(1.1 kpc) by 5 M_☉ pc⁻² decreases the expectation of Σ_>1.5(Z) by only 0.15 M_☉ pc⁻².

The derived profile of Σ_>1.5(Z) is shown in the lower panel of Figure 1. The same conclusions as before can be drawn: the results perfectly match the expectation of the visible mass alone, and the estimated local DM density is ρ_{☉, DM} = 0.0 ± 0.7 mM_☉ pc⁻³. Moreover, the discrepancy with the models comprising a DM halo is even more striking, since the OM model is 9σ higher than the derived solution, and the SHM and MIN models are excluded at the ∼7.5σ and ∼5σ level, respectively.

The calculation relies on a set of kinematical measurements, ten hypotheses, and three parameters. In this section, we will assess in more detail the reliability and robustness of the results, investigating if any of these input quantities can be the cause of the mismatch between observations and model expectations, and if the DM halo could pass undetected with our methods. Finally, alternative DM spatial distributions (non-spherical halo, dark disk, dark ring) will also be considered and compared to the observations.

4.1. The Incidence of the Parameters

A different definition of the solar Galactocentric distance has a negligible impact on the results. Both Σ(Z) and Σ_>1.5(4 kpc) decrease with R_☉, but they change by less than ±3.5 M_☉ pc⁻² and ±2.6 M_☉ pc⁻² (∼0.5σ at 4 kpc), respectively, when R_☉ is varied in the range 7.5–8.5 kpc. The precise value adopted for R_☉ is therefore irrelevant.

The effects of varying the thick disk scale height and length in the range 0.6–1.2 kpc and 3.0–4.6 kpc, respectively, are shown in Figure 2. Decreasing h_{Z, ρ} shifts the curve of Σ(Z) to higher values, but does not affect its slope noticeably. For example, Σ (1.5 kpc) agrees well with the expectation of the most massive DM halo model (OM) when assuming h_{Z, ρ} = 0.6 kpc, but at Z = 4 kpc they differ by ∼4σ, because the gradient of the observed curve is too shallow. In fact, even assuming this extremely low scale height, the less massive model (MIN) is still 2.5σ higher than the curve of Σ_>1.5(Z). On the contrary, varying the thick disk scale length affects the slope of Σ(Z), but does not add a lot of mass to the derived solution. Σ_>1.5(Z) increases with h_{R, ρ}, but even assuming h_{R, ρ} = 5 kpc, the MIN model is still 2σ higher than both Σ(Z) and Σ_>1.5(Z).

Zoom In Zoom Out Reset image size

While varying h_{Z, ρ} or h_{R, ρ} alone is not enough to reconcile the measurements with the presence of a classical DM halo, they can be changed simultaneously. Nevertheless, observational constraints prevent from obtaining a solution overlapping the expectations of DM halo models because even the less massive MIN model can be roughly matched by the observations only by assuming h_{Z, ρ} = 0.65 kpc and h_{R, ρ} = 4.7 kpc. Such a large-scale length was proposed by some authors (Ratnatunga & Freeman 1989; Beers & Sommer-Larsen 1995; Chiba & Beers 2000; Larsen & Humphreys 2003; Chang et al. 2011), but it was always associated with a higher scale height, as shown⁷ in Figure 3. The observations indicate that the required thin and extended thick disk is very unlikely. In conclusion, varying the three parameters involved in the calculations within the ranges allowed by the literature does not solve the problem of the missing DM in the volume under analysis.

Zoom In Zoom Out Reset image size

4.2. Non-flat Rotation Curve

The behavior of the Galactic rotation curve is still debated for R ⩾ 11 kpc (e.g., Binney & Dehnen 1997), but there is a general consensus that it is flat at the solar Galactocentric position. Indeed, recent observations find only a tiny negative gradient (Xue et al. 2008; Fuchs et al. 2009, $\partial \overline{V}/\partial R=-0.85$ and −0.006 ± 0.016 km s⁻¹ kpc⁻¹, respectively). Large deviations from a flat curve are therefore excluded, and theoretical models predict $\vert \partial \overline{V}/\partial R\vert \le 6$ km s⁻¹ kpc⁻¹ (e.g., Dehnen & Binney 1998; Olling & Merrifield 2000). Levine et al. (2006) also used a flat curve to measure the Galactic density distribution.

If the hypothesis (VIII) is dropped and a non-flat rotation curve is considered, the additional term

$$\begin{equation} \vspace*{-2pt}\frac{1}{\pi GR}\int _{0}^{Z}\frac{\partial \overline{V}}{\partial R}\overline{V} dz \end{equation} \tag{ 21 }$$

is added to the right-hand term of both Equations (13) and (14). In the following, we will assume V_LSR = 220 km s⁻¹, and the vertical shear of the thick disk rotational velocity $\overline{V}(Z)$ will be represented by the linear expression of Moni Bidin et al. (2012), that well approximates the underlying low-index power law in the Z-range under analysis. It is immediately evident from Equation (21) that a decreasing rotation curve at R_☉ (${\partial \overline{V}}/{\partial R}\le 0$ km s⁻¹ kpc⁻¹) introduces a negative term to the calculation, thus increasing the discrepancy between the observations and the expectations of the DM halo models. On the contrary, a higher Σ(Z) is obtained if the rotation velocity increases with R, and a certain amount of DM is thus allowed in the volume under analysis. However, ${\partial \overline{V}}/{\partial R}= 10$ km s⁻¹ kpc⁻¹ is required to match the minimum DM density deduced by the Galactic rotation curve (MIN model), and ${\partial \overline{V}}/{\partial R}= 16.5$ km s⁻¹ kpc⁻¹ to match the SHM. Such steep rotation curves are excluded by observations. Moreover, we are comparing the observations with models whose density was fixed to return a flat rotation curve, hence we step into the contradiction that, while assuming a steeper curve to increase the measured mass, we simultaneously increase the amount of DM required by the models. Among the possibilities offered by the literature, the solution closest to the expectations of a DM halo model is obtained assuming a rotation curve increasing its steepness from $({\partial \overline{V}}/{\partial R}){\vert _{Z=0}}=0$ to $({\partial \overline{V}}/{\partial R}){\vert _{Z=4 \,\mathrm{kpc}}}\approx 7$ km s⁻¹ kpc⁻¹, as modeled by Kalberla et al. (2007) for the Galactic gas. Even in this case, however, the MIN model is 2σ higher than both the measured Σ_>1.5(Z) and Σ (4 kpc), while the SHM is excluded at the 5σ level. In conclusion, assuming a non-flat rotation curve can alter the results, but insufficiently to justify the mismatch between the observations and the models.

4.3. Alternative Radial Profiles of the Dispersions

A popular model for the radial dependence of σ_U, alternative to the hypothesis (X), is the assumption of a constant Toomre Q-parameter. Its validity is controversial, and it is disfavored by the numerical integrations of Cuddeford & Amendt (1992). If the rotation curve is flat, the constant-Q model predicts that

$$\begin{equation} \sigma ^2_\mathrm{U}(R)\propto R^2 \exp \left(-\frac{2R}{h_{R,\rho }}\right) \end{equation} \tag{ 22 }$$

(Amendt & Cuddeford 1991), from which it immediately follows that

$$\begin{equation} \frac{\partial \sigma ^2_\mathrm{U}}{\partial R}=\sigma ^2_\mathrm{U}\cdot \left(\frac{2}{R}-\frac{2}{h_{R,\rho }}\right). \end{equation} \tag{ 23 }$$

When dropping the hypothesis (X) to assume the constant-Q model, a choice about the radial dependence of σ²_V and $\overline{UW}$ must also be done. However, they have a lower incidence on the results, and from Equation (23) it is easy to see that ∂σ²_U/∂R is almost indistinguishable in the two cases at the solar position, because R_☉ ≈ 2h_{R, ρ}. Hence, it makes little difference if they follow the trend of Equation (23), or a purely exponential decay. This is not the same for σ²_U because the second derivative ∂²σ²_U/∂R² enters into Equation (13).

Substituting the assumption (X) with Equation (23) for σ²_U leads to a solution identical to Equation (14), with k₁ replaced by

$$\begin{equation} k_1\mbox{$^\prime $}=\frac{15}{R_\odot \cdot h_{R,\rho }}-\frac{6}{R_\odot ^2}-\frac{6}{h_{R,\rho }^2}. \end{equation} \tag{ 24 }$$

The results are shown in Figure 4. The curve of Σ(Z) both increases its slope and shifts upward, and a certain amount of DM is allowed, but this is still insufficient when compared to the models: the expectation of the SHM is still ≈4.5σ higher than both Σ(Z) and Σ_>1.5(Z), while the discrepancy reduces to 2σ for the minimum-density MIN model. The assumption of the constant-Q model for the radial decay of σ²_U mitigates the gap between the observations and the models, but not enough to reconcile them.

Zoom In Zoom Out Reset image size

Neese & Yoss (1988) found that their measurements of σ_U throughout the Galactic disk could be fitted with a linear radial decay, with slope −3.8 ± 0.6 km s⁻¹ kpc⁻¹. The effects of dropping the assumption (X) in favor of linear radial profiles can be explored, although there is no theoretical support for this behavior. We will assume ∂σ_U/∂R from Neese & Yoss (1988), and ∂σ_V/∂R = −6.0 km s⁻¹ kpc⁻¹ from Casetti-Dinescu et al. (2011), whose results for σ_U are identical to those of Neese & Yoss (1988), although they do not claim that the underlying trend should be linear. Fixing ${\partial \overline{UW}}/{\partial R}$ is more problematic due to the absence of measurements in the literature, but from the last two terms of Equation (13) it is easy to see that a more negative value decreases the slope of Σ(Z). Therefore, we will adopt the most favorable case $\overline{UW} (R)$ = constant.

Assuming the radial linear trends given above, the calculation of Σ(Z) from Equation (13) is straightforward, and the results are shown in Figure 4. The curve of Σ(Z) has a much steeper slope than under the assumption (X), and Σ_>1.5(Z) matches the expectation of the MIN model, agreeing with the SHM model within ∼1σ. However, still Σ(Z) is offset by 10–15 M_☉ pc⁻² in the whole Z-range when compared to the MIN model, and 25–30 M_☉ pc⁻² with respect to the SHM curve. Thus, the theoretically unjustified linear radial profile of σ_U still requires an additional ad hoc correction, like a strong overestimate of the visible mass, or a thick disk scale height of h_{Z, ρ} = 0.7 kpc, and it is therefore unlikely.

4.4. A Flared Thick Disk

The flare of the Galactic thick disk is still an issue of debate, and the observations have not provided unique evidence of its existence. Siegel et al. (2002) and Du et al. (2006), for example, argue in favor of a flared thick disk with scale height increasing with R, while Cabrera-Lavers et al. (2007), supported by Bilir et al. (2008) and Yaz & Karaali (2009), proposed a flare of opposite sign. Given these contradictory results, it is most probably safe to assume in our calculations that the flare is negligible at the solar position, else a relevant variation of the thick disk scale height would have been detected beyond doubt.

If we assume a flared thick disk, thus dropping the assumption (IX), the scale height in Equation (9) becomes a function of R, and we have

$$\begin{equation} \frac{\partial \ln {(\rho)}}{\partial R}=-\frac{1}{h_{R,{\rm flare}}}= \left(-\frac{1}{h_{R,\rho }}+\frac{Z}{h^2_{Z,\rho }(R)}\cdot \frac{\partial h_{Z,\rho }}{\partial R}\right). \end{equation} \tag{ 25 }$$

The consequence of this new approach is the simple substitution of h_{R, ρ} with h_{R, flare} in Equation (13), plus the additional term

$$\begin{equation} -\int _{0}^{Z}\sigma ^2_U \cdot \frac{\partial }{\partial R}\left(\frac{1}{h_{R,{\rm flare}}}\right)dz. \end{equation} \tag{ 26 }$$

The new formulation, although analytically simple, hides many practical complexities: first, the expression for h_{Z, ρ}(R) must be accurate enough to provide a good approximation up to its second radial derivative. Additionally, the assumption (X) implicitly relies on the constancy of h_{Z, ρ}(R) (van der Kruit & Searle 1981, 1982), and in a flared disk the radial behavior of the dispersions is harder to model. Here we will therefore limit to the first-order approximation that the flare is locally small, hence well represented by a linear function of R (as proposed by Cabrera-Lavers et al. 2007), and not noticeably affecting the radial behavior of the dispersions. In this case, the additional term of Equation (26) is neglected, and a solution identical to Equation (14) is obtained, with k₁ and k₃, respectively, substituted by

$$\begin{eqnarray} &&k_1\mbox{$^\prime $}=\frac{1}{R_\odot \cdot h_{R,{\rm flare}}}-\frac{1}{h_{R,\rho }^2}+ \frac{1}{h_{R,\rho }}\cdot \Bigg(\frac{2}{R_\odot }-\frac{1}{h_{R,{\rm flare}}}\Bigg),\nonumber \\ \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&k_3\mbox{$^\prime $}=\frac{2}{h_{R,\rho }}+\frac{1}{h_{R,{\rm flare}}}-\frac{2}{R_\odot }. \end{eqnarray} \tag{ 28 }$$

The resulting profiles of the surface density are shown in Figure 5. The introduction of a flare in the formulation only changes the slope of the curve, which decreases in the case of a negative flare (∂h_{Z, ρ}/∂R < 0). A positive flare, however, does not return results that can be easily reconciled with the presence of a DM halo: Σ_>1.5(Z) can be forced to match the MIN model assuming ∂h_{Z, ρ}/∂R = 0.1, and the SHM model with ∂h_{Z, ρ}/∂R = 0.16, but even in this case Σ(Z) still remains 15–20 M_☉ pc⁻² lower than the expectations of this model in the whole z-range. As in the case of a linear radial decay of the dispersion (Section 4.3), an additional correction to the formulation is therefore required.

Zoom In Zoom Out Reset image size

4.5. Alternative Kinematical Results

We repeated the calculations replacing the kinematical measurements of Moni Bidin et al. (2012) with other works in the literature to check the incidence of the assumed kinematics on our results. Three previous investigations are suitable for our purposes: Fuchs et al. (2009), Casetti-Dinescu et al. (2011), and Smith et al. (2012). The results are shown in Figure 6.

Zoom In Zoom Out Reset image size

The data points of Fuchs et al. (2009) in the range Z = 0–0.8 kpc were linearly fitted to derive the vertical profiles of σ_U, σ_V, σ_W, and $\overline{UW}$ to be inserted in Equation (14). A quadratic fit led to negligible differences. We assumed h_{Z, ρ} = 0.3 kpc and h_{R, ρ} = 2.6 kpc (Jurić et al. 2008) because their sample mostly comprises thin disk stars. The curve of Σ(Z) thus derived is clearly unphysical. However, Binney (2009) has already shown that the kinematical data of Fuchs et al. (2009) lead to inconsistent results when used to constrain the mass distribution in the Galaxy. The most likely cause of the problem is that their sample is a mixture of old thin and thick disk stars, with the incidence of the thick disk increasing with Z, causing too steep a vertical gradient of the dispersions. On the contrary, Casetti-Dinescu et al. (2011) used a pure thick disk sample to measure the variation of the kinematics in the range Z = 0.3–2.2 kpc. Unfortunately, they did not give the Z-profile of $\overline{UW}$, but their measurement of the tilt angle α (strongly related to $\overline{UW}$) agrees with Moni Bidin et al. (2012). Hence we assumed the same profile of $\overline{UW} (Z)$ as before. In any case, the incidence of this term is very limited at lower Galactic heights. The use of the Casetti-Dinescu et al. (2011) data returns the same general results previously discussed: Σ(Z) well matches the expectations for the visible mass alone even in this case. Unfortunately, the results are weaker, both because of the larger observational errors, and because the expectations of the models are less distinct at lower Z. Thus, the less massive DM model is only 1σ higher than the calculated solution.

Smith et al. (2012) recently measured the change of the three-dimensional disk kinematics between Z = 0.5 and 1.8 kpc, in three metallicity bins. The metal-rich group ([Fe/H] ⩾−0.5) is probably a heterogeneous mix of thin disk sub-populations, and it is not suitable for our purposes. The metal-poor thick disk, dominating the sample with [Fe/H] ⩽−0.8, has a unique spatial distribution (Carollo et al. 2010) poorly studied in the literature. Hence, we adopted the results for the intermediate-metallicity stars (−0.8 ⩽ [Fe/H] ⩽−0.5), representative of a thick disk stellar population similar to that studied by Moni Bidin et al. (2012). To derive the vertical trend of the kinematical quantities, we linearly fitted their results excluding their nearest bin ($\overline{Z}$ = 0.7 kpc). We thus limited the incidence of a residual thin disk contamination, not efficiently excluded at low Z by the metallicity-only sample selection. Moreover, restricting the Z-range ensures that the linear fit is a good approximation of the underlying trend because close to the plane, as already discussed in Section 3, the vertical gradient of the dispersions is expected to vary rapidly with Z. We did not estimate an error on Σ(Z) in this case, because the uncertainties associated with the kinematical quantities are not well defined in our linear fit of three data points. As shown in Figure 6, the solution obtained with Smith et al. (2012) overlaps that derived by means of the results of Casetti-Dinescu et al. (2011), again roughly matching the expectations for the visible mass only. In conclusion, two more sets of independent kinematical results return results identical to ours when used in the calculation of Σ(Z), although with a lower significance.

4.6. DM Halo Detectability

The failure to detect the presence of a DM halo naturally raises the question whether this feature is detectable by our method. In fact, the calculation does not measure the mass directly, but a variation of the gravitational potential, from which the mass density is derived. Moni Bidin et al. (2010) noted that, if the DM halo was too uniform, extended, and poorly concentrated, it could cause a negligible change of the potential in the volume under analysis, thus resulting undetectable.

The gravitational potential can be separated at any point as the sum of the potential of the dark and visible matter, ϕ = ϕ_DM + ϕ_VIS. It is easy to see that the dynamical surface density calculated through Equation (1) can therefore be expressed as the sum of two contribution, Σ(Z) = Σ_{ϕ, DM}(Z) + Σ_{ϕ, VIS}(Z), where

$$\begin{eqnarray} &&\Sigma _{\phi,\mathrm{DM}}(Z)=\frac{1}{2\pi G} \left[\int _{0}^{Z} \frac{1}{R}\frac{\partial }{\partial R}\Bigg(R \frac{\partial \phi _\mathrm{DM}}{\partial R}\Bigg)dz +\frac{\partial \phi _\mathrm{DM}}{\partial Z}\right], \nonumber\\ \end{eqnarray} \tag{ 29 }$$

and analogously for Σ_{ϕ, VIS}(Z). The potential of a spherical NFW DM distribution can be written as

$$\begin{equation} \phi _\mathrm{DM}(s) = -g(c)\frac{GM_\mathrm{V}}{r_\mathrm{s} c}\frac{\ln {(1+cs)}}{s}, \end{equation} \tag{ 30 }$$

where r_s is a characteristic radial scale length, c is the concentration parameter (Lokas 2001), s = r/(r_sc) with $r=\sqrt{R^{2}+Z^{2}}$ the radial spherical coordinate, M_V is the total mass within the Viral radius r_V = r_sc, and g(c) = (ln (1 + c) − c/(1 + c))⁻¹. The two NFW models analyzed in this paper have M_V = 5.2 × 10¹¹ M_☉ and c = 16.8 (SHM model), and M_V = 6.85 × 10¹¹ M_☉ and c = 9 (N97 model). Inserting Equation (30) into Equation (29), we can estimate Σ_{ϕ, DM}(Z), i.e., the dynamical DM mass that can be detected by means of the Poisson equation between ±Z from the Galactic plane. The results of this calculation for the SHM and N97 models are shown in Figure 7. In both cases, the dynamical DM mass perfectly matches the DM mass expected by the models, indicating that all the mass of the DM halo enclosed between ±Z is detected by the use of Equation (1). It is therefore false that the small potential difference induced by the DM halo between 0 and 4 kpc from the plane can induce an underestimate of the physical mass in the volume under analysis.

Zoom In Zoom Out Reset image size

It can also be argued that, even if Equation (1) should not fail, the change of the stellar kinematics between Z = 1.5 and 4 kpc induced by the presence of a DM halo could be too small to be detected by the kinematical measurements adopted in our work. To check this possibility, we compared the vertical gradient of the potential of a DM halo, ∂ϕ_DM/∂Z, with the observed quantity:

$$\begin{equation} \frac{\partial \phi }{\partial Z}\Bigg {\vert }_\mathrm{obs}=\frac{\sigma ^2_W}{h_{Z,\rho }} -\frac{\partial \sigma ^2_W}{\partial Z} -\overline{UW}\cdot \Bigg{(}\frac{1}{R}-\frac{2}{h_{R,\rho }}\Bigg{)}. \end{equation} \tag{ 31 }$$

This is the vertical gradient of the gravitational potential estimated by means of the stellar kinematics through Equation (3). The comparison, relative to the SHM model, is shown in Figure 8. Clearly, the DM halo does not generate a strong vertical gradient of the potential: the surface density of the expected DM mass enclosed within ±4 kpc is similar to the measured one (∼60 M_☉ pc⁻²), while ∂ϕ_DM/∂Z is about a factor of two smaller than the observed potential gradient. Nevertheless, (∂ϕ/∂Z)|_obs is measured with an error of the order of 10% (∼200 km² s⁻² kpc⁻¹). This indicates that the presence of a DM halo would had not passed undetected because it would had caused a variation of the measured (∂ϕ/∂Z)|_obs about four to five times larger than the observational errors. In a forthcoming study, the orbit integration of a synthetic stellar population embedded in the Galactic potential will be used to test the reliability of our method, and to analyze in more detail the causes of the systematics that can arise.

Zoom In Zoom Out Reset image size

4.7. Dark Halo Flattening

So far, we have considered only spherical models for the DM halo, but the expected amount of DM in the volume under analysis changes noticeably if its shape is modified. The effects of varying the flattening of the models SHM and MIN are shown in Figure 9. The models were calculated from Equations (19) and (20), where R_c, ρ_c, and ρ_{☉, DM} were multiplied by 1/q to keep the resulting Galactic rotation curve the same (Olling 1995; Olling & Merrifield 2001). Clearly, oblate models (q < 1) imply a larger quantity of DM in the solar vicinity, and they depart further from the observations. The expected surface density, on the other hand, decreases at increasing q, and highly prolate models approach the observational results. A perfect match requires a zero-density halo and is therefore never reached for any finite q, but a lower limit can be derived. We thus conclude that strongly prolate models are required, because within 2σ (95% confidence level) we have q ⩾ 2 for the MIN model, and this lower limit is much higher for more common, higher-density models, for example, q ≳ 4 for the SHM.

Zoom In Zoom Out Reset image size

4.8. Other Dark Matter Structures

We have compared the observations with common models of DM halos, but other spatial distributions have been presented in the literature, like a DM disk and a ring (e.g., Kalberla et al. 2007). Nevertheless, these features were usually proposed in addition and not as an alternative to a spheroidal halo, thus enhancing the expected local DM density.

The presence of a DM ring at R ≫ R_☉ cannot be detected with our observations, but most of the mass of a DM disk would be included in the volume under analysis, and thus detectable by our method. The existence of such a feature, first proposed by Lake (1989), is of particular interest, because it is a natural expectation of current ΛCDM models of Galactic formation (Purcell et al. 2009). In Figure 10 our results are compared to the expectations of the massive DM disk model of Kalberla (2003), the thick, low-density model of Purcell et al. (2009), and the thinner and denser ones from Read et al. (2008). The comparison shows that the dark disk proposed as an alternative to a DM halo (Kalberla 2003) is excluded by the observations with very high significance (9σ) and, as discussed in the previous sections, this mismatch cannot be corrected by simply altering one parameter or one assumption. On the contrary, the local density of the DM disk models proposed by Read et al. (2008) and Purcell et al. (2009) is low, and they cannot be ruled out. Some constraints can be derived, though. The curves of both Σ(Z) and Σ_>1.5(Z) indicate that, within 2σ, the thick DM disk of Purcell et al. (2009, scale height ∼4.6 kpc) must have a local density ρ_{☉, DM} < 3 mM_☉ pc⁻³, while the upper limit for the thinner Read et al. (2008) model (scale height of 2.4 kpc) is ρ_{☉, DM} < 5 mM_☉ pc⁻³. However, it must be taken into consideration that these results are sensitive to a change of a parameter or an assumption as discussed in Sections 4.1–4.4. In any case, these DM disk models do not sustain the Galactic rotation curve, and their presence in addition to a DM halo would enhance the discrepancy between the observations and the expected DM density in the volume under study.

Zoom In Zoom Out Reset image size

The measurement of the mass surface density at the solar Galactocentric position between 1.5 and 4 kpc from the Galactic plane accounts for the visible mass only. The DM density in the solar neighborhood, extrapolated from the observed curve of Σ(Z), is ρ_{☉, DM} = 0 ± 1 mM_☉ pc⁻³, at variance with the general consensus that it must be in the range 5–13 mM_☉ pc⁻³ (e.g., Weber & de Boer 2010; Garbari et al. 2011). Our recent measurements of the thick disk kinematics were used in the calculation, but the observed lack of DM is independent of this choice, because very similar results can be obtained by means of other kinematical results in the literature. The calculation relies on three input parameters and ten assumptions, but the observations cannot be reconciled with the DM halo modes modifying one of them. Altering them at will introduces enough freedom to force the solution to match the expectations of the most preferred model, but in this case an exotic series of unlikely hypotheses must be invoked. For example, a very thin thick disk (h_{Z, ρ} = 0.7 kpc) either very extended in the radial direction (h_{R, ρ} = 4.6 kpc) or strongly flared at the solar position (∂h_{Z, ρ}/∂R = 0.1) makes the solution coincide with the minimum DM local density deduced from the Galactic rotation curve, but the observational constraints exclude or disfavor these scenarios. On the other hand, the expected visible mass strikingly matches the observations without any effort, by use of the most probable assumptions. This coincidence lends weight to the interpretation of these results because it is easy to obtain an unphysical solution if one or more wrong hypotheses are being made. For example, Moni Bidin et al. (2010) showed that the results obtained under the assumption of a cross-term $\overline{UW}$ symmetric with respect to the Galactic plane, as an alternative to our hypothesis (IV), violate two minimum requirements: the surface density must at least account for the known visible matter, and it cannot decrease with Z. Thus, the excellent agreement between the measured mass and the visible mass is unlikely to have been obtained by pure chance. We also demonstrated that our method cannot fail to detect the presence of a classical DM halo because it causes a noticeable change in the stellar kinematics, one order of magnitude larger than the observational errors.

The only viable solution to reconcile the observations with the models is the assumption of a highly prolate DM halo, that can sustain a flat rotation curve with a negligible density at the solar position. Observational constraints on the DM halo shape are scarce and often controversial and, while spherical structures are usually preferred (e.g., Ibata et al. 2001; Majewski et al. 2003; Johnston et al. 2005; Fellhauer et al. 2006), a prolate spheroid (q = 5/3) was invoked by Helmi (2004) and Law et al. (2005). Our observations suggest that q ⩾ 2 is required even in the lower bound case of the least massive model, for agreement within 2σ. However, very prolate structures are atypical in cold DM simulations, which have problems in reproducing them (e.g., Dubinski & Carlberg 1991). Thus, it must still be proven that a DM halo with a high flattening parameter is fully compatible with the current ΛCDM paradigm.

A dark spheroidal component is required to sustain the Galactic rotation curve, observationally confirmed to be flat for R ⩾ 5 kpc (Kalberla et al. 2007; Sofue et al. 2008; Xue et al. 2008). In the presence of the visible mass only, Newtonian dynamics would predict a steep Keplerian fall-off. Very noticeably, the calculation returns the unphysical result of a total surface density decreasing with Z if a rapidly decreasing curve is assumed. Hence, a flat rotation curve is required, while failing to detect the DM necessary to sustain it. This apparent contradiction is actually a confirmation that the calculation is reliable because it is consistent with an observationally proven fact, although not with its expected explanation.

In conclusion, the observations point to a lack of Galactic DM at the solar position, contrary to the expectations of all the current models of Galactic mass distribution. A DM distribution very different to what it is today accepted, such as a highly prolate DM halo, is required to reconcile the results with the DM paradigm. The interpretation of these results is thus not straightforward. We believe that they require further investigation and analysis, both on the observational and the theoretical side, to solve the problems they present. We feel that any attempt to further interpret and explain our results, beyond that presented in this paper, would be highly speculative at this stage. Future surveys, such as Gaia, will likely be crucial to move beyond this point. However, as our results currently stand, we stress that, while numerous experiments seek to directly detect the elusive DM particles, our observations suggest that their density may be negligible in the solar neighborhood. This conclusion does not depend on the cause of this lack of DM at the solar position. For example, if our results are interpreted as evidence of a highly prolate cold DM halo with q ⩾ 2, this would have a local density lower than 2 mM_☉ pc⁻³, i.e., more than a factor of four lower than what usually assumed in the interpretation of the results of these experiments.

It is clear that the local surface density measured in our work, extrapolated to the rest of the Galaxy, cannot retain the Sun in a circular orbit at a speed of ∼220 km s⁻¹. A deep missing mass problem is therefore evidenced by our observations. Indeed, we believe that our results do not solve any problem, but pose important new ones.

The authors warmly thank the referee C. Flynn for his extensive work of revision of our results, and for his detailed criticisms. C.M.B. and R.A.M. acknowledge support from the Chilean Centro de Astrofísica FONDAP No. 15010003, and the Chilean Centro de Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL PFB/06. R.S. was financed through a combination of GEMINI-CONICYT fund 32080008 and a COMITE MIXTO grant. C.M.B. also thanks F. Mauro for useful discussions. All authors acknowledge partial support from the Yale University/Universidad de Chile collaboration. The SPM3 catalog was funded in part by grants from the US National Science Foundation, the Yale University, and the Universidad Nacional de San Juan, Argentina. We warmly thank W. F. van Altena, V. I. Korchagin, T. M. Girard, and D. I. Casetti-Dinescu for their help and suggestions.

Facilities: Du Pont (ECHELLE) - Du Pont Telescope, Magellan:Clay (MIKE) - Magellan II Landon Clay Telescope, Euler1.2m (CORALIE) - The Swiss Euler 1.2m Telescope, Max Planck:2.2m (FEROS) - Max Planck Gesellschaft 2.2 meter Telescope at European Southern Observatory