A new determination of the local dark matter density from the kinematics of K dwarfs

We apply a new method to determine the local disc matter and dark halo matter density to kinematic and position data for \sim2000 K dwarf stars taken from the literature. Our method assumes only that the disc is locally in dynamical equilibrium, and that the 'tilt' term in the Jeans equations is small up to \sim1 kpc above the plane. We present a new calculation of the photometric distances to the K dwarf stars, and use a Monte Carlo Markov Chain to marginalise over uncertainties in both the baryonic mass distribution, and the velocity and distance errors for each individual star. We perform a series of tests to demonstrate that our results are insensitive to plausible systematic errors in our distance calibration, and we show that our method recovers the correct answer from a dynamically evolved N-body simulation of the Milky Way. We find a local dark matter density of {\rho}dm = 0.025+0.014-0.013 M\odotpc^{-3} (0.95+0.53-0.49 GeV cm^{-3}) at 90% confidence assuming no correction for the non-flatness of the local rotation curve, and {\rho}dm = 0.022+0.015-0.013 M\odotpc^-3 (0.85+0.57-0.50 GeV cm^{-3}) if the correction is included. Our 90% lower bound on {\rho}dm is larger than the canonical value typically assumed in the literature, and is at mild tension with extrapolations from the rotation curve that assume a spherical halo. Our result can be explained by a larger normalisation for the local Milky Way rotation curve, an oblate dark matter halo, a local disc of dark matter, or some combination of these.


INTRODUCTION
The local dark matter density is an average over a small volume, typically a few hundred parsecs, around the Sun. It provides constraints on the local halo shape and allows us to predict the flux of dark matter particles in laboratory detectors. The latter is required to extract information about the nature of a dark matter particle from such experiments, at least in the limit of a few tens to hundreds of detections (Peter 2011). The Galactic halo shape can be constrained by combining two methods of determining the local dark matter density. Firstly, one can infer it from the Galactic rotation curve (ρ ext dm ). This requires an assumption about the shape of the Galactic halo (typically spherical; e. g. Sofue et al. 2009;Weber & de Boer 2010;Catena & Ullio 2010). Secondly, one can calculate the dark matter density locally from the vertical kinematics of stars near the Sun (ρ dm ) (e. g. Bahcall e-mail: silvia@physik.uzh.ch 1984b; Holmberg & Flynn 2000). If ρ dm < ρ ext dm , this suggests a prolate dark matter halo for the Milky Way; while ρ dm > ρ ext dm , could imply either an oblate halo or a dark disc (Lake 1989; Read et al. 2008;Read et al. 2009).
Determining the local matter density from the kinematics of stars in the Solar Neighbourhood has a long history dating back to Oort (1932Oort ( , 1960 in the 1930's. Oort used the classical method of solving the combined Poisson-Boltzmann equations for a sample of stars, assumed to be stationary in the total matter distribution of the disc. He found 50% more mass than the sum of known components. A more modern study by Bahcall (1984b) introduced a new method that described the visible matter as a sum of isothermal components. He also found dynamically significant dark matter in the disc 1 (Bahcall 1984a). Using faint K dwarfs at the South 1 We should be careful about what we mean by 'dark matter in the disc'. Early studies like Oort (1932) were typically interested in missing disc-like matter (a 'thin dark disc'); more mod-arXiv:1206.0015v2 [astro-ph.GA] 23 Jul 2012 Galactic Pole, Bahcall et al. (1992) confirmed his earlier result that more than 50% of the mass was dark, although with a lower statistical significance. However, the early studies by Oort (1932Oort ( , 1960 and Bahcall (1984b,a) assumed that different tracers could be simply averaged to form a single tracer population. Kuijken & Gilmore (1989a) demonstrated that the two samples of F stars analysed by Bahcall (1984b) were not compatible (i.e. they had different spatial density distribution, but no evidence for a difference in their kinematics) and therefore should not be averaged. They reanalysed the K giant sample used by Bahcall (1984a), assigning more realistic errors to the density profile and using a more detailed fit to the velocity data, finding a value of total matter density compatible with the observed one. They concluded that the determination of the local volume density remained limited by systematic and random errors with the available data.
With the launch of the ESA satellite Hipparcos (1997), the kinematics and position of tracer stars were measured with much higher accuracy. The improved distance measures give a much more accurate measurement of the local luminosity function, so that the total amount of visible matter can be better estimated as well. The latest dynamical measurements of the local density of matter -ρtot -from Hipparcos data show no compelling evidence for a significant amount of dark matter in the disc (Creze et al. 1998;Holmberg & Flynn 2000). Holmberg & Flynn (2000) found ρtot = 0.102 ± 0.01M pc −3 , with a contribution of about 0.095M pc −3 in visible matter, consistent with the 's value.
In addition to the local volume density, several authors have calculated the local surface density of gravitating matter, probing up to larger heights above the disc plane (typically ∼ 1 kpc; e.g. Kuijken & Gilmore 1989b;, 1991Holmberg & Flynn 2004). Using faint K dwarfs at the South Galactic Pole, and using a prior from the rotation curve, , 1991 find ρ KG dm = 0.010±0.005 M pc −3 , consistent with that expected from the rotation curve assuming a spherical Galactic dark matter halo 2 (e. g. Sofue et al. 2009;Weber & de Boer 2010;Catena & Ullio 2010). A similar result was found in the post-Hipparcos era by Holmberg & Flynn (2004). Recently, Moni Bidin et al. (2012) have estimated the surface density using tracers at heights 1.5 < z < 4 kpc above the disc, making a rather stronger claim (incompatible with the earlier results of Kuijken & Gilmore (1991) and Holmberg & Flynn (2004)) that there is no dark matter near the Sun. However, Bovy & Tremaine (2012) demonstrate that this result is erroneous and owes to one of ten assumptions used by Moni Bidin et al. (2012) being false. Furthermore, Sanders (2012) estimate that the velocity dispersion gradients derived by ern studies try to constrain a significantly more extended dark matter halo that has a near-constant dark matter density up to ∼ 1 kpc. Even the 'dark disc' predicted by recent cosmological simulations (Read et al. 2008;Read et al. 2009) is sufficiently hot that its dark matter distribution is approximately constant up to ∼ 1 kpc. Throughout this paper when we talk about 'dark matter in the disc' we refer to a constant density dark matter component within the disc volume. 2 Note that this consistency with the rotation curve is somewhat circular since this is input as a prior in their analysis.
Moni Bidin et al. (2012) could be biased by up to a factor of two, which would also significantly alter their determination of ρ dm .
With next generation surveys round the corner (e.g. Gaia; Jordan 2008), a significant improvement in the number of precision astrometric, photometric and spectroscopic measurements is expected. For this reason, Garbari et al. (2011) (hereafter Paper I) revisited the systematic errors in determining ρ dm from Solar Neighbourhood stars; these will likely soon become the dominant source of error, if they are not already. We were the first to use a high resolution Nbody simulation of an isolated Milky Way-like galaxy to generate mock data. We used these mock data to study a popular class of mass modelling methods in the literature that fit an assumed distribution function to a set of stellar tracers (Holmberg & Flynn 2000;Binney & Tremaine 2008). We found that realistic mixing of stars due to the formation of a bar and spiral arms (similar to those observed in the Milky Way) breaks the usual assumption that the distribution function is separable, leading to systematic bias in the recovery of ρ dm . We then introduced a new method that avoids this assumption by fitting instead moments of the distribution function (i.e. that solves the Jeans-Poisson equations). Our Minimal Assumption method (or MA method) uses a Monte Carlo Markov Chain technique (hereafter MCMC) to marginalise over remaining model and measurement uncertainties. Given sufficiently good data, we showed that our method recovers the correct local dark matter density even in the face of disc inhomogeneities, non-isothermal tracers and a non-separable distribution function.
In this article, we apply our MA method to real data from the literature. The key advantages of our new method over previous works are that: (i) we use a 'minimal' set of assumptions; (ii) we use a MCMC to marginalise over both model and measurement uncertainties; and (iii) we require no prior from the Milky Way rotation curve as has been commonly used in previous works , 1991Holmberg & Flynn 2000, 2004. This latter means that we can compare our determination to that derived from the rotation curve to constrain the Milky Way halo shape. Our method requires at least one equilibrium stellar tracer population with known density fall off ν(z) and vertical velocity dispersion σ 2 z (z), both as a function of height z. The requirements for a suitable sample of stellar tracers are that: (i) they are in dynamical equilibrium with the Galactic potential (i.e. they must be sufficiently dynamically old to have completed many vertical oscillations through the Galactic plane); (ii) they are available in sufficient numbers to give good statistical precision; (iii) they have reliable distances and vertical velocities vz; (iv) the sample completeness needs to be sufficiently well understood in order to measure the density fall off as a function of the distance z from the disc plane; and (v) they extend up to 2-3 times the disc scale height (in order to break a degeneracy between the disc and dark matter densities; see Paper I). While full six dimensional phase space information is now available for a large number of stars (e.g. RAVE Steinmetz 2003;Steinmetz et al. 2006;Zwitter et al. 2008;and SEGUE Yanny et al. 2009), these surveys are magnitude rather than volume complete, with additional survey selection effects based on colour. This makes it difficult to reliably estimate ν(z) for a given tracer population. For this reason, we return to the volume com-plete K-dwarf data from  for our disc tracers -the 'KG' data. These data consist of a photometric sample of 2016 K dwarf stars, complete in the z-range ∼ 0.2 − 1.5 kpc, with a spectroscopic sample of 580 K dwarfs (most of which are included in the photometric catalogue). We use data from Hipparcos and SEGUE (Kotoneva et al. 2002;Zhang 2012) to perform a new photometric distance measurement for each K dwarf star. We model the local gravitational potential using the baryonic mass distribution of the Galactic disc by Flynn et al. (2006). This article is organised as follows. In Section 2, we present the K dwarf data from  (hereafter KG989II) and describe our new distance determinations (2.1). In Section 3, we summarise our MA mass modelling method (3.1) and, for comparison, the method adopted by KG89II (3.2). In Section 3.3, we test both methods on a mock data set derived from a dynamically evolved N-body simulation. In Section 4, we apply our MA method to the KG data and present our results. Finally in Section 5, we summarise and present our conclusions.

DATA
KG89II present a catalogue of 2016 K stars with photometry in B and V bands, and another of 580 K dwarfs (most of which are also included in the photometric catalogue) including radial velocities, at South Galactic pole (figure 1). The stellar density fall off of these tracers was derived from star counts. At large z, the mean metallicity of the stars is known to decrease below the Solar Neighbourhood value. Such a gradient translates into an absolute magnitude gradient, since the position of the main sequence in the colourmagnitude diagram changes with the metallicity: metal poor dwarfs are fainter than metal rich ones at the same temperature or colour (the opposite is true for giant stars). So, if there is a vertical metallicity gradient, the photometric parallaxes used for the derivation of the density fall off -ν(z) -will be systematically wrong as one moves away from the plane. Unfortunately, the metallicity could not be measured directly for these stars, so KG89II derived the density fall off using an assumed constant metallicity gradient for the K dwarfs. They considered two different gradients to estimate the magnitude of the uncertainties, namely d[Fe/H]/dz = 0 (constant metallicity) and d[Fe/H]/dz = −0.3 dex kpc −1 (in their analysis KG89II consider this latter as the fiducial metallicity gradient for K dwarfs).
The tracers' vertical distance determination is fundamental for our analysis. Twenty years after KG89II's study, we can re-calibrate the distances for these stars using modern survey data to estimate the metallicity distribution function of K dwarfs at different z, and Hipparcos parallaxes to calibrate the photometric distances. Our distance recalibration procedure is described, next.

A new distance determination for the K dwarfs
To calculate the vertical distances z for KG89II's sample, we must derive a relationship between the metallicity [Fe/H], the vertical distance z and the absolute magnitude MV of K dwarf stars in the disc. Since the metallicity is not included in KG89II's catalogue, we can only hope to derive a distance distribution function P * (z) for each star of the sample, based on an observed metallicity distribution function for these stars. We consider two different catalogues of K dwarfs with distances and metallicity for our calibration. The first catalogue by Kotoneva et al. (2002) consists of 431 K dwarfs from Hipparcos, representing a complete catalogue of the metal content in nearby K dwarfs extending up to z ∼ 100 pc. The vertical distances z for these stars are very accurately determined from Hipparcos parallaxes. The second catalogue by Zhang (2012) contains 5000 SEGUE K dwarfs spanning a much wider range of z, namely between 300 and 2000 pc. However, in this case, the distance determination for these stars is much less certain: the distance errors are about 10%.
We combine these two catalogues to build the metallicity distribution function (MDF) Q([Fe/H](z), z) for K dwarfs. Comparing the SEGUE metallicity distribution function at z = 500 pc with the MDF from Kotoneva et al. (2002) for z < 100 pc (see the red and black solid histograms in figure 2, respectively), we notice that the two MDFs have very similar shape, but shifted means. There is a vertical metallicity gradient from 0 to 500 pc of ∼ −0.4 dex kpc −1 (corresponding to a shift of −0.2 dex; see figure 2 black dotted histogram). At larger height than this, the gradient is weaker: the SEGUE MDF at 1 kpc (dashed red histogram) is similar to the one at 500 pc. We adopt the Kotoneva et al. (2002) MDF in the Galactic plane, then we apply a linear metallicity gradient of ∼ −0.4 dex kpc −1 between 100 pc and 500 pc to match the SEGUE MDF at z = 500 pc; we extend this shifted Kotoneva et al. (2002)'s MDF up to z = 800 pc, and we adopt the SEGUE MDF for z > 800 pc, as shown in figure 3. We explore an alternative Q([Fe/H](z), z) for the K dwarfs in Appendix A.
With the metallicity distribution function    Kotoneva et al. (2002) and SEGUE (Zhang 2012 The best-fit parameters are [−5.795, 27.92, 0.1291 where V is the apparent magnitude and AV is the extinction; we use the value AV = 0.062 mag from Schlegel et al. (1998), given the mean Galactic coordinates of KG89II's data. The vertical distance z for a single star is then: where l and b are the Galactic longitude and latitude. We know l, b and B − V for each star of KG89II's sample, so the only free parameter is [Fe/H]. This means that the vertical distance for each star will be given by a probability distribution P * (z) corresponding to a metallicity distribution P * ([Fe/H]) for that star: In practice, this equation must be solved iteratively because [Fe/H] is itself a function of z through equation 3. The iterative process proceeds as follows: (i) We start the first iteration by assuming that the distance distribution function of a single star is P * (z) = 1 for all the possible z([Fe/H]) calculated using equation 3.
(ii) The MDF marginalised over z for a single star is given by: where Q([Fe/H](z), z) is the observed MDF at each z for all K dwarfs as described previously.
(iii) A new P * (z) is calculated through equation 4, using the P * ([Fe/H]) just computed.
(iv) We restart steps (ii) to (iii) to calculate P * ([Fe/H]) with the new P * (z) until it converges. For all stars, we obtained convergence in less than 5 iterations.
The density fall off of the photometric sample and the velocity dispersion of the spectroscopic one, obtained with our new distance estimates, are shown in figure 5. For our analysis, we use the density and the velocity dispersion profiles only over the range 200 < z < 1200 pc (red dashed lines). This assures that our sample is volume complete and avoids significant contamination by K giant stars . The corresponding quantities computed by KG89II, assuming a metallicity gradient of −0.3 dex/kpc, are plotted as a comparison (empty circles). Our new density profile and velocity dispersion do not differ greatly from those of KG89II; they are compatible within the quoted errors. In Appendix A, we explore the effect of a rather extreme variation in the assumed MDF -ignoring Kotoneva    (2002) and SEGUE data -finding that our results are not sensitive to plausible changes in our distance calibration. The referee of this article pointed out that the study of Kotoneva et al. (2002) has been updated by Casagrande et al. (2007). The two studies are very much compatible, but the scatter in equation 1 of 0.03 mag becomes 0.27 mag when the newer data are used. We tested the impact of this larger scatter in magnitude on the distance calibration, finding that the density and velocity dispersion profile remain unchanged in the range of z of interest, with a negligible increase in the uncertainties (see Appendix A).

The MA method
The MA method presented in Paper I uses the Poisson-Jeans system to predict the density fall off of a tracer population in a given gravitational potential. The comparison between this predicted density fall off and the observed one allows us to constrain the gravitational potential and, consequently, the underlying dark matter distribution.
Here we summarise the basic equations; for a detailed description of the MA method see Section 2.1 of Paper I.
The MA method is based on three main assumptions: (i) The system is in equilibrium (steady state assumption).
(ii) The dark matter density is constant over the range of |z| considered.
(iii) The 'tilt' term 1 R ∂ ∂R Rνσ 2 Rz in the cylindrical Jeans equation: is negligible compared to all other terms. Here ν, σ 2 z and σ 2 . Upper panel: K dwarf stellar density profile (filled circles with error bars) derived from a Monte Carlo sampling of P * (z) for each star. As a comparison, the density profile (assuming a metallicity gradient of −0.3 dex/kpc) from KG89II is plotted as empty circles. Lower panel: The similarly derived vertical velocity dispersion as a function of z (filled circles with error bars). The corresponding determination by KG89II is represented by the empty circles. In both panels, the two red dashed lines show the range of z considered in our analysis, over which the photometric sample is volume complete and we avoid significant contamination from K giants stars.
are the number density and the velocity dispersion components of a tracer population moving in potential Φ.
With these assumptions, the Jeans equation becomes a function only of z and we can neglect the other two Jeans equations in R and θ: Solving this equation for a single tracer population, we obtain its density ν(z) at each height z: Given the density at the midplane ρs,j(0) and the vertical velocity dispersion σ 2 z,j (z) as a function of z for each of the gas and stellar populations in the local disc, we can model the full disc density distribution as a superposition of such elements: In Paper I, we showed that accurate measurement of the vertical velocity dispersion of the tracers σ 2 z (z) is crucial, however in the mass modelling we can assume that all the visible matter components are isothermal -i.e. σ 2 z,j = σ 2 z,j (0) -and equation 9 simplifies to: The Poisson equation then determines the potential Φ from the density. In cylindrical polar coordinates, this is given by: with: where ρ dm (R) is the halo mass density (following assumption (ii), this is assumed to be independent of z in the volume considered); and Vc(R) = (R∂Φ/∂R) 1/2 is the (total) circular velocity at a distance R (in the plane) from the centre of the Galaxy. For a flat rotation curve, the second term vanishes and ρ eff dm (R) = ρ dm (R). The rotation curve correction can be calculated from the Oort constants A and B (Binney & Merrifield 1998): We solve equations 10 and 11 numerically for a given tracer population, adopting the following procedure: (i) We make initial trial guesses for ρs,j(0), ρ dm , and the vertical velocity dispersion for the visible matter component in the plane σ 2 z,j (0). (ii) We solve equation 10 to obtain Φ(z) and its first derivative ∂Φ ∂z , with Φ(0) = ∂Φ ∂z 0 = 0. (iii) We insert this result into equation 8 for the vertical density fall off of the tracers νp(z) and we compare this with the observed distribution ν(z) to obtain a goodness of fit.
This procedure requires many input parameters, such as the normalisations and dispersions of each of the disc components, and the vertical dispersion profile of the tracers (typically poorly constrained). To explore the parameter space and marginalise over the uncertainties, we use a Monte Carlo Markov Chain (MCMC) method.
In Paper I we showed that our method is able to recover the correct value of the local dark matter density, even in presence of large visible matter density fluctuations due to the spiral arms. Notice that the MA requires no prior from the Milky Way rotation curve, as has been commonly used in previous works; this means that we can compare our determination to that derived from the rotation curve to constrain the Milky Way halo shape.

Application of the MA method to real data
When we apply the MA method to real data, we must deal with distance and velocity uncertainties, and account for survey geometry and/or the sample completeness. In particular, the K dwarf data from KG89II are not assigned distances, but instead, as described in Section 2.1, a z probability distribution function P * (z) for each star of the sample. In addition, the vertical velocities for each star are measured with an uncertainty of about 0.5 − 1 km/s. To marginalise over these uncertainties, we proceed in the following way: (i) For each model n of the MCMC, we select a different vertical distance z * n for each star in the sample, according to its z distribution function P * (z). We make sure that for each star included in both the spectroscopic and the photometric sample we pick a unique z * n value. For each star of the spectroscopic sample we also draw a vertical velocity value v * z,n from a Gaussian distribution, according to its velocity error bar.
(ii) We bin the data in z to construct the observed density fall off νn(z) and velocity dispersion σ 2 z,n (z) for the tracers, selecting stars with z between 0.2 and 1.2 kpc. The velocity dispersion is calculated using the velocity scale algorithm described in Beers et al. (1990).
(iii) We use the trial guesses for ρs,j,n(0), σ 2 z,j,n (0) and ρ dm,n to solve equation 10 and 11 to obtain Φn(z) and its first derivative ∂Φn ∂z . (iv) We insert this result and the vertical velocity dispersion of the tracers σ 2 z,n (z) into equation 8 to predict the tracers' density fall off. When applying the MA method to the real data, we add the visible matter surface density Σs as a further constraint. For each model of the MCMC we compute Σs as (v) We compare the predicted density fall off νn,p(z) with νn(z), the predicted visible matter surface density Σ p s,n with the observed one (Σs) and we calculate the corresponding χ 2 n , accepting or rejecting the model n.
(vi) When a model is accepted, we restart from (i) with the following model (n + 1) and so on, exploring the whole parameter space.

The KG method
KG89II used the K dwarf data to calculate the total surface density at the Sun position Σ(R ). Their approach (the KG method) is similar to the method adopted later by Holmberg & Flynn (2000, 2004 (the HF method), that we analysed in detail in Paper I.
Instead of measuring the vertical velocity dispersion of the tracers as a function of z to predict the density fall off, the HF method uses the tracers' velocity distribution function in the mid-plane of the Galactic disc fz(vz,0); assuming that the distribution function is separable, Holmberg & Flynn (2000, 2004 integrate this distribution function over z-velocities to predict the density fall off of the tracers (see Section 2.2 of Paper I for more details): where fz(z, vz) = fz(Ez) and Ez = 1 2 v 2 z +Φ(z) is the vertical energy. This equation can be written as: where vz,0 is the vertical velocity at the Galactic mid-plane (z = 0). The MA demands only that the tilt term in the Jeans equation (6) is small with respect to the other terms, the HF method requires the stronger assumption that the z-motion (and so the distribution function) is completely separable from the motion in the radial and azimuthal directions; this latter implies that the tilt term is exactly zero.
The HF approach has the advantage of exploiting the whole available information about the shape of the velocity distribution function of the tracers. However, we demonstrated in Paper I that when the separability of the distribution function is not fulfilled, the HF method leads to biased results. Using a high resolution simulation of a Milky Way like galaxy, we showed that the onset of spiral arms and a bar can cause significant radial mixing that breaks the separability of the motion in the z and R directions, violating this key assumption of the HF method. This effect becomes increasingly important with height z. It is important to notice that the HF method can not be corrected for this bias since the separability of the potential (and of the distribution function) lies at the heart of the method. By contrast, in the MA method, the separability of the potential enters only in the neglected tilt term of the Jeans equation (that is assumed to be small as compared to the other terms). If the tilt term is for some reason large -i.e. the radial derivative of the density weighted tilt of the velocity ellipsoid is large -a correction can straightforwardly be applied to our MA method. However, we expect the tilt term to be small (see Paper I and Binney & Tremaine 2008), and so we do not consider such a correction in this paper.
KG89II's approach relies on the same key assumption about the distribution function as the HF method. In most studies, the density ν(z) of the tracers is known to better precision than the velocity distribution fz(vz, z). For this reason, KG89II work in the opposite direction with respect to Holmberg & Flynn (2000, 2004 and predict fz(vz, z) from the observed ν(z). Applying an inverse Abel transform to equation 15, they obtain: so there is a unique relation between ν(Φ) and fz(Ez). Notice that fz(Ez) depends on ν(Φ(z)) only at large z, where the potential exceeds Ez, i.e. beyond z = Φ −1 (Ez). Thus an additional key advantage of the KG method is that one can model the potential at large distances from the Galactic plane, ignoring the detailed distribution of matter at small z. KG89II parameterised the gravitational potential Φ(z) above the bulk of the disc matter (where it is sensitive only to the total surface density of gravitating matter) as: where D is the disc scale height, K is proportional to the total disc surface density Σ(R ), and F ∝ ρ eff dm (the effective halo density). KG89II used a range of Galactic mass models (calculated using different values of the disc mass M , the radial disc scale-length R d , the circular velocity Vc(R ) and Sun position R ) to ensure consistency with the Galactic rotation curve (assuming a spherical Milky Way halo) and therefore to obtain a relation between F and K. Note that already this is different from our MA approach where we use no information about the rotation curve to constrain our mass models.
Given the observed space density of a tracer population ν(z) and a set of gravitational potential models Φ(z), one can solve equation 17. To reduce the noise in the differential of ν(z), KG89II fitted it with a double exponential. KG89II then used the derived fz(Ez) for each potential model Φ(z) to compute the likelihood of the spectroscopic sample: where the product is over all stars in the spectroscopic sample, and select the potential parameters that maximise this likelihood function L.
The KG method, like the HF method, uses the full shape of the observed velocity distribution function, maximising the use of the available information. It is also convenient because it does not require a detailed model of the gravitational potential or an accurately measured tracer density fall off close to the Galactic plane. However, its drawback is that, like the HF method, it relies on a key assumption that the vertical distribution function is only a function of Ez. In the following section we test, using the high resolution N-body simulation described in Paper I, how this assumption affects the result derived using the KG method.

Testing the methods using an N-body simulation
Before applying the MA method to the real K dwarf data, we use the most dynamically evolved stage of the simulation described in Paper I as a mock data set, to test the effect of the velocity errors and of the asymmetric distance distribution function P * (z) on the MA method's result. Then we use the same mock data to probe how the non-separability of the tracers' distribution function affects the KG method.
The mock data consist of a high resolution (30 million ρ dm Figure 6. Results for recovering ρ dm and ρs from the simulation. Top left panel: The recovered dark matter density (red filled circles) for all 8 volumes analysed around the disc (90% and 68% confidence intervals are shown). The actual value of the dark matter density is marked as a blue filled circle. Top right panel: The density of models explored by the MCMC projected onto ρ dm -ρs space for the 90 • volume. The blue dot shows the true value for ρ dm and ρs; the diagonal cyan blue line shows the total matter density; and the red dot shows the median recovered ρ dm and ρs, with 90% errors marked. Bottom left panel: Histogram of the recovered ρ dm from the MCMC mode ensemble for this same volume. The striped grey area is the 90% confidence interval. The cyan dashed line and shaded area give the actual value of ρ dm with error bars. Bottom right panel: Testing the effect of assuming a separable distribution function. Each dot in this plot shows the likelihood Ln of each MCMC model calculated assuming a separable distribution function (see equation 17; the plot shows − log(Ln), so the more likely models have a lower ordinate value). The red dashed line shows a fit to the points. Notice that assuming that the distribution function is separable produces a bias towards low ρ dm (the true distribution function for the simulation is not separable; see figure 7). This bias effect was observed in all the volumes considered.
disc star particles and 15 million halo dark matter particles) N-body simulation of an isolated Milky Way like galaxy. The initial conditions were built to contain some thousand stars in the volume size required for our analysis. With the dynamical evolution of the simulation, the disc developed a bar and spiral arms. For more details about the simulation features and how it compares to the real Milky Way, see Section 3.1 of Paper I.
We consider several volumes in the disc of the simulated galaxy at a distance of R = 8.5 kpc from the centre. The MA method solves the Jeans equation for a onedimentional slab (equation 7), so the radial size of the volumes, ∆R = 0.25 kpc, is chosen to fulfil this approximation, but still contain an enough large number of stars (see Section 3.3.1.1 of Paper I for more details).
As described in Section 3.1.1, we assign a different velocity value v * z,n and a different z * n to each star at every iteration n of the MCMC. In the application of the MA method to the simulation, the velocity values are drawn from a Gaussian distribution centred on the true velocity value and with a width of 1 km/s, while the z values are selected from a lognormal distribution around the true value. Because of the numerical resolution of the simulation, we cannot fit the density profile up to 1.2 kpc, since at such height we quickly run out of star particles and the velocity dispersion is poorly measured, so we use stars with 0.2 z 0.75 kpc. For the simulation data, we model the visible mass distribu-tion as a single component, characterised by its mass density ρs,j = ρs(0) and its velocity dispersion on the midplane σ 2 z,j = σ 2 z (0). We let the dark matter density freely vary between 0 and 0.2 M pc −3 , and the other parameters -ρs(0) and σ 2 z (0) -vary inside their error bars. We adopted Poisson errors for the velocity dispersion σ 2 z (0) and the current uncertainties on the total visible matter density in the plane (i.e. 0.014 M pc −3 , see Section 4.1) for ρs(0). We include the rotation curve correction term, which can be easily computed for the simulation, in the calculation. We test convergence of our MCMC chains by starting with ρ dm seeded at two different values (namely ρ dm = 0 and ρ dm = 0.2 M pc −3 ) and running until the two chains are statistically indistinguishable.
In figure 6, the results for the MA method are shown. The upper left panel shows the results for eight different volumes around the simulated disc. Notice that in all cases, the mean correct answer (blue filled circles) is recovered within the 90% confidence interval, while for four out of eight of the patches (with a fifth at 45 • extremely close) it is recovered with the 68% confidence interval. This is consistent with our confidence intervals having the meaning of a purely statistical error, despite each patch being systematically different (each patch samples a different region of the disc with different local dynamics). The remaining three panels focus on the results for the volume at 90 • . The top right panel shows the density of models explored by the MCMC projected onto ρ dm -ρs space; the bottom left panel shows a histogram of the dark matter density for all the models explored by the MCMC; and the bottom right panel explores the effect of assuming a separable distribution function (of which more, next).
In Paper I, we noticed that, at this evolved stage of the simulation, the distribution function f (vR, v θ , vz, R, z) of the tracers is not separable. This leads to the HF method producing biased results -either underestimating or overestimating the local dark matter density (see Section 3 of Paper I for more details).
The separability of the distribution function lies at the heart of the KG method too. To reproduce a KGlike method, we consider all νn(z) and model potentials Φn(z) explored by our MA method MCMC chain. For each model in the chain, we then calculate a distribution function fz,n(Ez(vz, Φ)) through equation 17 and use it to compute the likelihood Ln of the velocity data via equation 19.
In practice, the integral in the denominator of equation 19 is calculated numerically from E min z = 1 to E max z = 7000 km 2 s −2 , which is chosen to avoid the divergence at Ez = 0 and to ensure we cover all energies of interest (the contribution of the high energy tail of fz,n(Ez) is negligible with respect to the low energy part; see figure 7).
In the bottom right panel of figure 6, the likelihoods Ln of the velocity data are plotted against the corresponding values of the local dark matter density for the different MCMC models ρ dm,n . This panel shows that there is an anti-correlation between the computed likelihood Ln and the corresponding value of the local dark matter density ρ dm,n : the likelihood of the velocity data is larger (i.e. − log(Ln) is lower) for gravitational potential models with low ρ dm,n ; this means that we expect the KG method to artificially favour low dark matter density values. For all the explored volumes around the disc, we always obtain this same anti- The distribution function -f -of the tracers in the simulation. The black line shows distribution function measured directly from the simulation averaged over a volume 8.375 < R < 8.625 kpc; 0.2 < z < 0.75 kpc. The orange dashed line shows fz(Ez) predicted from the density fall off ν(z) via equation 17, assuming the correct value of ρs and ρ dm in the computation of the gravitational potential. The blue and red dashed lines show fz(Ez) predicted from ν(z), assuming ρ dm = 0 and twice ρ dm , respectively. All distribution functions are normalised by the integral of fz(Ez) between the minimum and the maximum Ez of the star particles in the volume considered.
correlation, so the bias on ρ dm will have always the same sign.
To understand this effect, in figure 7 we show the distribution function fz(Ez) calculated from the same density profile ν(z), but using three potentials with different values of ρ dm , namely ρ dm = 0 (blue dashed line), the true value of ρ dm (orange dashed line) and twice the true ρ dm (red dashed line); the black line represents the actual distribution function of the stars in the volume. The potential corresponding to the true value of the dark matter density does not predict the distribution function correctly, while with ρ dm = 0 we obtain a better agreement with the measured fz(Ez). It is clear from this plot that the likelihood of models with low dark matter density, calculated through equation 19, will be higher than the true model. Therefore the KG method will be biased towards low ρ dm .

Measuring the local matter and dark matter density
In the previous section, we showed that the MA method is able to recover the correct dark matter density within our 90% confidence interval, even in presence of asymmetric distance errors, velocity uncertainties and a non-separable distribution function. We now apply the MA method to the real K dwarf data, proceeding as described in Section 3.1.1. The mass distribution in the Galactic disc is modelled as a superposition of 15 isothermal components, listed in Table 1. As parameters to fit in the MCMC, we use the local dark matter density ρ dm , the total visible density in the midplane ρs(0), and the relative fractions of the visible components ρs,j(0)/ρs(0) and their velocity dispersions in the midplane σz,j. We allow the densities and the velocity dispersions of the different components to vary within their measured uncertainties (the errors for each component are given in Table 1). We let the total visible density in the plane ρs(0) vary within its observed range: ρs(0) = 0.0914 ± 0.0140 M pc −3 (extrapolated from where: and where the sum is extended to all the bins and ∆νn,i are the uncertainties on the density fall off. The rotation curve correction term can be calculated from the Oort constants A and B. To determine the Oort constants, we must use stellar tracers that are well-mixed. The most recent estimates of A and B from F giants (Branham 2010) and K-M giants (Mignard 2000) from Hipparcos give A = 14.85 ± 7.47 km s −1 kpc −1 and B = −10.85±6.83 km s −1 kpc −1 and A = 14.5±1.0 km s −1 kpc −1 and B = −11.5 ± 1.0 km s −1 kpc −1 , respectively. Averaging these two values we obtain a correction term of −0.0033 ± 0.0050 M pc −3 . We test for convergence of the MCMC by starting several chains at different initial values of all the parameters and running until they are statistically equivalent (after removing an initial burn-in phase of 100 accepted models for each chain).
The results for the MA method applied to the real K dwarf data are shown in figure 8. The upper panel shows the density of the models explored by the MCMC (grey contours) in the ρs − ρ dm plane; the median with 90% errors is shown by the black dot and corresponds to ρ dm = 0.025 +0.014 −0.013 M pc −3 (0.95 +0.53 −0.49 GeV cm −3 ) 3 ; adding the rotation curve correction, we obtain ρ dm = 0.022 +0.015 −0.013 M pc −3 (0.85 +0.57 −0.50 GeV cm −3 , see the red dot in figure 8).
We expect the MA method to primarily constrain the total matter density in the plane ρs + ρ dm , which means that we should see an oblique degeneracy between ρs and  Table 1. The disc mass model taken from Flynn et al. 2006. For each component in the table, we give the local mass density in the midplane ρ(0) in M pc −3 ; the total column density Σ in M pc −2 ; and the vertical velocity dispersion σ z,j (0) in km s −1 . Uncertainties on the densities are assumed to be 50% for all the gas components (indicated with * ) and 20% for all of the stellar components. The largest uncertainties come from the gas that remains poorly constrained (compare, for example, compilations in Flynn et al. 2006, ? andFerrière (2001)). For the thick disc, the column density is well known, while the velocity dispersion and the volume density are poorly known such that they should have larger error bars. However, these two quantities are essentially nuisance parameters for our analysis here. Since they anticorrelate and -as pointed out by Kuijken & Gilmore (1989b) the local gravitational potential is mainly constrained by the column density, we simply assume small errors for both here such that the integrated column agrees with the observed value. ρ dm in the ρs − ρ dm plane, as was observed for the simulation (see figure 6). However, unlike the simulation that has only one visible matter component in the disc, our real-data mass model comprises some 15 separate components with different scale heights. This introduces new freedoms that wash out the oblique degeneracy in the ρs − ρ dm plane. If we plot instead, however, the total surface density of visible matter in the disc (bottom panel of figure 8) the degeneracy is once again clearly visible as a diagonal elongation of the MCMC model density contours. Notice that many of the models explored by the MCMC lie outside of the range given by the Flynn et al. 2006 mass model (grey striped band). This simply means that the prior we placed on Σs is not very strong. However, the full area explored is in good agreement with the more conservative measurements of the total visible matter surface density by Kuijken & Gilmore (1991) (hereafter KG91), namely Σs = 48 ± 8 M pc −2 ; and by Flynn & Fuchs (1994), Σs = 49 ± 9 M pc −2 . Even if we include only models that lie within the grey striped region, our recovered ρ dm is little affected.
In Appendix B, we test the robustness of our result, exploring how it changes if one considers either only the low z bins (0.2 < z < 0.7 kpc) or only the high z bins (0.6 < z < 1.2 kpc). We find that the low z data do not provide any information about ρ dm , but they still favour low Σs compared to the prior we imposed. The low z bins present a noisier and not monotonically increasing velocity dispersion, however they are better constrained and dominate the χ 2 fit. Using the high z data, we lose information about the disc and Σs settles into the centre of its prior distribution; this leads to a systematically lower ρ dm , since the sum of the two is well constrained. This tells us that the origin of our high ρ dm is the slightly lower Σs required by the velocity dispersion data near the plane.
We also tested the effect of changing the assumed errors on the baryonic mass model. We first reduce them to an optimistic 10% error for the stellar normalisations and 30% for the gas normalisations; this has little effect on the resulting determination of ρ dm because sufficient freedom remains in our mass model to allow degeneracies between ρ dm and the baryonic components (this reflects the broken degeneracy seen in figure 8). We then increased the errors by removing the prior on ρs(0) altogether. This also has a small effect on our recovered ρ dm . This confirms the results shown in Appendix B that it is the velocity dispersion data, not our prior that are constraining our mass model. The low z ( 500 pc) data constrain the surface density profile of the disc, while the high z data ( 500 pc) are dynamically sensitive to dark matter.
The green horizontal stripe in the upper panel of figure  8 marks the value of ρ dm we extrapolate from KG89II and KG91. KG91, using the same K dwarf data analysed in this paper, determine the total dynamical surface density up to 1.1 kpc: Σ dyn = 71 ± 6 M pc −2 . If we subtract from this the contribution of the observed visible matter Σs = 48 ± 8 M pc −2 , we can calculate ρ dm , assumed to be constant in the range 0 < z < 1.2 kpc, as: This gives: ρ KG dm = 0.010 ± 0.005 M pc −3 . Our new result from our MA method is in tension with ρ KG dm , obtained from the same data set. This could owe either to our different distance calibration or to our new MA method that does not require any assumption about the separability of the distribution function (see Section 3.3). In figure 5, we already showed that our new distance calibration does not significantly affect the velocity dispersion and the density fall off of the K dwarfs. In addition, in Appendix A, we explore the effect of using a constant metallicity gradient for the K-dwarfs of −0.3 dex kpc −1 , exactly as assumed by KG89II. This metallicity distribution is not compatible with modern data; we only use it to illustrate the sensitivity of our results to the MDF of the K dwarf stars, and to fully understand why our determination of ρ dm is larger than that of KG91. Using the KG89II's MDF, our recovered value of ρ dm is slightly smaller and therefore in better agreement with KG89II and KG91. However, our median recovered value remains significantly larger than the upper bound of the KG91 result. This suggests that our new distance determinations are not the primary reason for the systematic shift.
In our tests on the N-body simulation (Section 3.3), we showed that, when the distribution function of the tracers is not separable, the method adopted by KG89II leads to a systematic underestimate of ρ dm . In the lower panel of figure 9, we plot the likelihood Ln of each model explored by our MCMC -calculated through equation 19 -against the corresponding value of ρ dm . Unlike the similar plot for our simulation data (figure 6, bottom panels), there is now a significant vertical dispersion in the models. This owes to the increased freedom present in our 15-parameter mass model for the real-data. However, the highest likelihood models (the bottom envelope of points in the plot) show a similar trend as seen for the simulation data: higher likelihood models have systematically smaller ρ dm . We conclude that the primary difference for our larger value of ρ dm as compared to KG89II is that our MA method requires no assumption about the separability of the distribution function.
In the upper panel of figure 9, we plot a histogram of ρ dm from all the models explored by the MCMC. The striped area is the 90% confidence interval (corresponding to the black dot of figure 8); the result including the rotation curve correction is shown by the red error bar. Notice that our 90% lower bound is larger than the Standard Halo Model 4 typically assumed in the literature (marked by the vertical dashed orange line). For a comparison, we plot the ranges of ρ dm at the solar radius obtained by Iocco et al. (2011), combining microlensing and rotation curve measurements, and using different halo models: the blue error bar corresponds to a spherical halo, while the cyan and purple bars correspond to oblates halos with potential flattening q = 0.9 and q = 0.7, respectively. The magenta bar represents the dark matter density in presence of a dark disc, contributing 0.25 − 1.5 times the dark matter (spherical) halo density, as predicted by Read et al. (2009).
From figure 9, we can see that our recovered density is in mild tension with the result for a spherical Milky Way halo. Moving to an oblate halo significantly reduces this tension, however a flattening of q = 0.7 is likely inconsistent with measurements of the halo shape from the Sagittarius stream of stars (e.g. Ibata et al. 2001). If we wish to explain our median value for ρ dm that is very much larger than the canonical SHM value assumed in the literature, we require a local disc of dark matter that raises ρ dm without significantly altering the rotation curve. Interestingly, our median value is in excellent agreement with the range of dark discs predicted for our Galaxy by Read et al. (2008); Read et al. (2009).

DISCUSSION AND CONCLUSIONS
We have presented a new measurement of the local matter and dark matter densities from the kinematics of K dwarf stars near the Sun. We presented a new photometric distance calibration for the the K dwarf data of KG89II (the KG data), derived using modern survey catalogues and the Hipparcos satellite data. We then used these data as tracers of the local gravitational potential to calculate the visible (ρs) and dark matter (ρ dm ) densities at the solar position R and the surface density of the Milky Way disc up to 1.1 kpc above the plane (Σs).
To determine ρ dm and ρs, we applied our new mass modelling method (presented already in Paper I) that relies on a minimum set of assumptions (the MA method) to . Upper panel: A histogram of the recovered ρ dm from our MCMC chains for the MA method applied to the real K dwarf data. The striped grey area is the 90% confidence interval. The orange dashed line is the SHM value of ρ dm ; the blue, cyan and purple error bars correspond to the value of ρ dm obtained by Iocco et al. (2011) from a combination of microlensing and rotation curve data, with a spherical halo (potential flattening q = 1) and two oblate halos, with q = 0.9 and q = 0.7, respectively. The magenta error bar is the value of ρ dm expected if the Milky Way has a dark disc contributing 0.25 − 1.5 times the density of the (spherical) halo. The red error bar corresponds to our result after adding the rotation curve correction. Lower panel: The effect of assuming a separable distribution function. Each dot shows the likelihood of a given MCMC model calculated assuming a separable distribution function. Notice that the assumption of separability biases the result towards low ρ dm . The orange and the green dashed lines have the same meaning as in the upper panel of 8.
the rejuvenated KG data. The key advantages of our new method are that: (i) we do not require any hypothesis about the shape of the tracers' velocity distribution function; (ii) we use a MCMC to marginalise over uncertainties in the distances and velocities of the tracer stars, and the underlying baryonic mass model for the visible disc; and (iii) we require no prior from the Milky Way rotation curve as has been commonly used in previous works. This latter means that we can compare our determination to that derived from the rotation curve to constrain the Milky Way halo shape. We used a dynamically evolved high resolution N-body simulation of a Milky Way-like galaxy as a mock data set to test our MA method, finding that we could correctly recover ρ dm and ρs within our 90% confidence interval (for eight sample Solar neighbourhood-like volumes) even in the face of disc inhomogeneities, non-isothermal tracers, asymmetric distance errors and a non-separable tracer distribution function. Furthermore, we confirmed the result from our Paper I that assuming a separable distribution function (as has been typically done in the modern literature) leads to a biased determination of ρ dm .
Applying our MA method to the K dwarf data, we obtain a new measurement of the local dark matter density: ρ dm = 0.025 +0.014 −0.013 M pc −3 (0.95 +0.53 −0.49 GeV cm −3 ); which, adding a correction for the local non-flatness of the rotation curve correction term ( −0.0033 ± 0.0050, see Sections 3.1 and 4), gives: ρ dm = 0.022 +0.015 −0.013 M pc −3 (0.85 +0.57 −0.50 GeV cm −3 ). Our new value is systematically larger than the results from KG89II and KG91 derived from the same data. We show that this primarily owes to our new MA modelling method (and the fact that it does not assume a separable distribution function for the tracers); our new distance determination for the K dwarfs plays a more minor role. At the same time we determine a value of the local visible matter density of ρs = 0.098 +0.006 −0.014 M pc −3 that largely reflects the prior from our baryonic mass model.
Our error bars are larger than is often quoted in the literature, however they reflect the full combination of model systematic, measurement and statistical uncertainties. Other recent determinations either rely on the rotation curve and therefore a strong assumption about the Milky Way halo shape, or require a large number of assumptions with associated (and typically unmodelled) systematic errors.
In addition to measuring ρ dm and ρs, we also obtain an estimate of the baryonic disc mass up to z = 1.1 kpc above the disc plane: Σs = 45.5 +5.6 −5.9 M pc −2 at 90% confidence. This is slightly lower than the mean of the prior from our baryonic mass model: Σvis = 49.4 ± 4.6 M pc −2 . Splitting the number into the contribution from stars and stellar remnants: Σ * = 33.4 +5.5 −5.2 M pc −2 and gas: Σg = 12.00 +1.9 −2.0 M pc −2 , we see that our model favours slightly lower surface density in both the gas and the stars than the mean of our priors (Σ obs g = 13.3 ± 3.4 M pc −2 , Σ obs * = 36.1 ± 3.0 M pc −2 ).
It is this tendency for our models to favour lower disc surface density that leads to our high median value for ρ dm (see figure 8 and Appendix B). Unfortunately, current estimates of the stellar and gaseous inventory in the Solar neighbourhood are too uncertain to confirm or rule out our favoured Σs (e.g. Bovy et al. 2012;Ferrière 2001).
Our median value of the local dark matter density is larger at 90% confidence than the Standard Halo Model value of ρ SHM dm = 0.008 M pc −3 (0.30 GeV cm −3 ) usually adopted in the literature. This could be a statistical fluctuation; in one out of eight patches in our simulated mock data, our method overestimated ρ dm by ∼ 90%. However, if our high median value is confirmed by future data then it has some interesting implications. Firstly, it is particularly important for direct detection experiments because it implies a larger flux of dark matter particles and therefore a greater chance of detection. Secondly, our result is at mild tension with the value of ρ ext dm extrapolated from the rotation curve measurements, assuming a spherical dark matter halo. This suggests that the halo of our Galaxy is oblate and/or that we have a disc of dark matter, as predicted by recent cosmological simulations (see upper panel of figure 9).

APPENDIX A: TESTING THE ROBUSTNESS OF THE MDF
The most uncertain quantity in our re-analysis of KG89II's data is the variation of the K dwarfs' metallicity distribution function with z, Q([Fe/H](z), z). In this appendix, we investigate how the adopted metallicity distribution function affects the result of our analysis by exploring a different model for this function. We use the gradient adopted by KG89II, i.e. −0.3 dex kpc −1 , and set the mean metallicity to 0 at z = 0 (see upper left panel of figure A1). This MDF is not compatible with modern metallicity data for the K dwarfs (see figure 2); we use it simply to illustrate our sensitivity to the assumed MDF, and to aid comparisons with the earlier KG89II results.
In figure A1, we show the velocity dispersion (upper right panel) and the density fall off (lower left panel) of the ttracers derived using the above MDF. In the lower right panel of figure A1, the recovered dark and visible matter density are shown. The value of ρ dm obtained is slightly lower than that derived using our default MDF, namely ρ dm = 0.022 +0.013 −0.012 M pc −3 , or ρ dm = 0.018 +0.014 −0.013 M pc −3 including the rotation curve correction. However, our median value for ρ dm is still high and the overall result remains in tension with the SHM value (even when using this in-correct MDF). We conclude that our results are robust to plausible variations in the assumed MDF.
In addition to the above, we also tested the impact of a different extinction value used in equation 2 on our results. We chose AV = 0 and AV = 0.1. For AV = 0 (AV = 0.1) the distances are slightly overestimated (underestimated) with respect to the extinction values considered in this article. This does not affect much the density fall off, but it translates mainly in a slightly flatter (steeper) velocity dispersion. This leads to a small decrease (increase) of the recovered value of ρ dm . However the impact of the extinction is very small compared to the previous test shown in this Appendix.
Finally we tested the impact on our distance calibration of the larger scatter in equation 1 obtained using the more modern K dwarf catalogue from Casagrande et al. (2007) instead of Kotoneva et al. (2002). The two studies are very much compatible, but the scatter in equation 1 of 0.03 mag increases to 0.27 mag when the newer data are considered. We used the 104 K dwarf stars from Casagrande et al. (2007) to build the relationship between MV and (B − V , [Fe/H]), similar to equation 1, but using the polynomial down to a2,1(B − V ) 2 [Fe/H]. Unlike Kotoneva et al. (2002), the new catalogue also provides the error on the parallax, therefore we use a χ 2 instead of a simple linear-least-square fit. To account for the 0.27 mag uncertainties on MV , we convolved the distance probability distribution function with a Gaussian of width σ = 13% (corresponding to 0.27 mag in MV ). The resulting density and velocity dispersion profiles are unchanged in the range of z of interest, with only a negligible increase in the uncertainties. Unlike in the case of the extinction test, where distances were (slightly) systematically overestimated or underestimated, the magnitude scatter has only the effect of slightly increasing the random errors. The effect is weak so as long as a sufficiently large number of stars is used per bin, and the scatter is not too large. For this reason our determination of the local dark matter density is not affected by the increased scatter found by Casagrande et al. (2007). SHM median (90%) median (90%) + RC KG Figure A1. Left upper panel: The MDF adopted by KG89II: a constant metallicity gradient with height of 0.3 dex kpc −1 normalised to a metallicity of 0. dex at z = 0 kpc. Right upper panel: Stellar density profile derived from Monte Carlo sampling the probability distribution of z (dots with error bars). As a comparison, the density profile from KG89II is plotted as empty circles. The two red dashed lines show the completeness range. Lower left panel: The vertical velocity dispersion as a function of z. The black dots with error bar are derived from the probability distribution of z. Lower right panel: The projection of MCMC models onto the ρ dm -ρs plane for our MA method applied to the K dwarf data, using the MDF shown in the top left panel to determine the distances. The colours and the symbols are as in figure 8.

APPENDIX B: EXPLORING THE ROBUSTNESS OF OUR ρ dm DETERMINATION
In this Appendix, we explore the robustness of our determination of ρ dm by analysing a low z (0.2 < z < 0.7 kpc) and a high z (0.6 < z < 1.2 kpc) subset of the KG data. The results are shown in Figure B1. If we consider only the low z data (left panel), there is no information about the dark matter density. However, the data still favour low Σs compared to our prior. If only the high z data are used (right panel), we lose the information about the disc and Σs settles into more or less the centre of its prior distribution. This leads to a systematically lower ρ dm . The above suggests that the origin of our high median ρ dm is the lower Σs favoured by the velocity dispersion data close to the plane. 0.2 < z < 0.7 kpc 0.2 < z < 1.2 kpc 0.6 < z <  Figure B1. The recovered dark matter density as a function of the visible matter surface density Σs. Left panel: considering only low z bins (0.2 < z < 0.7 kpc); central panel: considering the full z range: 0.2 < z < 1.2 kpc (as lower panel of figure 8); right panel: considering only high z bins (0.6 < z < 1.2 kpc). The meaning of the symbols is the same as in figure 8. The striped grey area is the range in the visible matter surface density determined by Flynn et al. (2006), that we used as a (weak) prior.