Issue 
A&A
Volume 571, November 2014



Article Number  A92  
Number of page(s)  13  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201424478  
Published online  14 November 2014 
Weighing the local dark matter with RAVE red clump stars
^{1}
Observatoire astronomique de Strasbourg, Université de Strasbourg, CNRS,
UMR 7550, 11 rue de l’Université,
67000
Strasbourg,
France
email:
olivier.bienayme@unistra.fr
^{2}
Mount Stromlo Observatory, RSAA, Australian National
University, Weston Creek,
ACT
2611
Canberra,
Australia
^{3} Institute for Computational Astrophysics, Dept of Astronomy
& Physics, Saint Marys University, Halifax, NS, BH3 3C3, Canada
^{4}
Jeremiah Horrocks Institute, University of Central
Lancashire, Preston,
PR1 2HE,
UK
^{5}
Institute of Astronomy, Cambridge University,
Madingley Road, Cambridge
CB3 0HA,
UK
^{6}
Astronomisches RechenInstitut, Zentrum für Astronomie der
Universität Heidelberg, Mönchhofstr. 12–14, 69120
Heidelberg,
Germany
^{7}
Sydney Institute for Astronomy, School of Physics A28, University
of Sydney, Sydney
NSW
2006,
Australia
^{8} INAF Osservatorio Astronomico di Padova, 36012 Asiago ( VI),
Italy
^{9} Department of Physics and Astronomy. University of Victoria,
Victoria, BC. Canada V8P 5C2
^{10}
Department of Physics, Macquarie University,
Sydney
NSW
2109,
Australia
^{11}
Research Centre for Astronomy, Astrophysics and Astrophotonics,
Macquarie University, Sydney, NSW
2109,
Australia
^{12}
Mullard Space Science Laboratory, University College
London, Holmbury St
Mary, Dorking,
RH5 6NT,
UK
^{13}
Department of Physics and Astronomy, Padova
University, Vicolo dellOsservatorio
2, 35122
Padova,
Italy
^{14}
LeibnizInstitut für Astrophysik Potsdam (AIP),
An der Sternwarte 16,
14482
Potsdam,
Germany
^{15}
Australian Astronomical Observatory, PO Box 915, North Ryde
NSW
1670,
Australia
^{16}
Johns Hopkins University, Homewood Campus, 3400 N Charles Street,
Baltimore, MD
21218,
USA
^{17}
Faculty of Mathematics and Physics, University of
Ljubljana, 1000
Ljubljana,
Slovenia
Received: 26 June 2014
Accepted: 17 September 2014
We determine the Galactic potential in the solar neigbourhood from RAVE observations. We select red clump stars for which accurate distances, radial velocities, and metallicities have been measured. Combined with data from the 2MASS and UCAC catalogues, we build a sample of ~4600 red clump stars within a cylinder of 500 pc radius oriented in the direction of the South Galactic Pole, in the range of 200 pc to 2000 pc distances. We deduce the vertical force and the total mass density distribution up to 2 kpc away from the Galactic plane by fitting a distribution function depending explicitly on three isolating integrals of the motion in a separable potential locally representing the Galactic one with four free parameters. Because of the deep extension of our sample, we can determine nearly independently the dark matter mass density and the baryonic disc surface mass density. We find (i) at 1 kpc K_{z}/ (2πG) = 68.5 ± 1.0 M_{⊙} pc^{2}; and (ii) at 2 kpc K_{z}/ (2πG) = 96.9 ± 2.2 M_{⊙} pc^{2}. Assuming the solar Galactic radius at R_{0} = 8.5 kpc, we deduce the local dark matter density ρ_{DM}(z = 0) = 0.0143 ± 0.0011 M_{⊙}pc^{3} = 0.542 ± 0.042 Gev cm^{3} and the baryonic surface mass density Σ_{bar} = 44.4 ± 4.1 M_{⊙}pc^{2}. Our results are in agreement with previously published K_{z} determinations up to 1 kpc, while the extension to 2 kpc shows some evidence for an unexpectedly large amount of dark matter. A flattening of the dark halo of order 0.8 can produce such a high local density in combination with a circular velocity of 240 km s^{1}. It could also be consistent with a spherical cored dark matter profile whose density does not drop sharply with radius. Another explanation, allowing for a lower circular velocity, could be the presence of a secondary dark component, a very thick disc resulting either from the deposit of dark matter from the accretion of multiple small dwarf galaxies, or from the presence of an effective “phantom” thick disc in the context of effective galacticscale modifications of gravity.
Key words: galaxies: kinematics and dynamics
© ESO, 2014
1. Introduction
The complexity of the structure, dynamics, and history of our Galaxy is progressively unveiled with the recent advent of numerous large surveys. The access to positions, velocities, and chemical abundances with reasonable accuracy for large samples of stars allows us to explore the detailed properties of our own Galaxy. In the near future, covering a huge Galactic volume with an unprecedented accuracy, Gaia observations will revolutionize our understanding of Galactic formation and evolution (Perryman et al. 2001).
For the time being, we concentrate on the question of the dynamical estimate of the mass distribution in the solar neighbourhood, the K_{z} problem (where K_{z} is the vertical force perpendicular to the Galactic plane). We note that a major difficulty, in any survey analysis, is the identification of systematic errors when they are an order of magnitude below the dispersions of the measurement errors. By selecting a sample of red clump (RC) stars from the RAVE survey (RAdial Velocity Experiment survey, Kordopatis et al. 2013), we can drastically improve the measurements needed for the K_{z} problem in terms of number or precision. In particular, for the first time, we succeed in separating the force contributions from the Galactic discs and from the dark halo, and we find that the local density of dark matter is higher than previously expected.
Besides a better description of the local morphology of our Galaxy, our determination of a large local density of dark matter has implications on the morphology of the dark matter component(s) and also has implications for the terrestrial experiments of direct detection of the dark matter. The local density of dark matter can also provide a test for the MOND effective theory (Milgrom 1983; Famaey & McGaugh 2012).
Our new dynamical determination of the Galactic potential and K_{z} force perpendicular to the Galactic plane is based on RAVE observations (DR4, Kordopatis et al. 2013) towards the South Galactic Pole (SGP). This paper is an extension of the previous works published by Soubiran et al. (2003), Siebert et al. (2003), Bienaymé et al. (2006), and Soubiran et al. (2008), which probed the properties of RC stars within 100 pc of the Sun and at larger distances towards the North Galactic Pole (NGP).
The advantage of RC stars is that their distances are quite accurately deduced from photometric measurements as long as such stars have been identified to belong to the red clump, usually from their colour and gravity (see Cannon & Lloyd 1969, for an early recognition of the red clump).
From previous RAVE data releases (Steinmetz et al. 2006; Zwitter et al. 2008; Siebert et al. 2011b), RC stars were already used to probe the Galactic structure, for instance for the identification of stellar populations (Veltz et al. 2008), for a first measurement of the velocity ellipsoid tilt (Siebert et al. 2008), or as a probe of the local 3D velocity field in a large volume of the solar neighbourhood (Siebert et al. 2011a; Williams et al. 2013).
Here, we consider a sample of about 4600 RC stars, mainly selected in the direction of the SGP up to distances of 2 kpc from the Galactic plane.
A novelty of this work, setting aside the size of the sample of stars with accurate distances, is our ability to measure the vertical potential further away from the plane, up to a vertical distance of 2 kpc. From this, we deduce the local surface mass density and we constrain the shape of the total vertical mass distribution.
In this paper, the methods developed for the data analysis and modelling are the standard methods, except for the introduction of a separable Stäckel potential with an explicit third integral of the motion to allow for the analysis at z distances higher than 1 kpc. In addition, the selection function of RAVE observations of RC stars is well defined and accurately determined. The completeness of observed RAVE RC stars towards the SGP is 83% at 700 pc, 66% at 1.5 kpc, and 20% at 2 kpc.
The history of the K_{z} measurements covers decades of Galactic astronomy, and various summaries can be found in previous publications (Read 2014). We just mention here that a turning point was the advent of the Hipparcos satellite observations, which were decisive by increasing the amount of accurate data and critical by allowing us precise calibrations of the absolute magnitude of stars. In particular, this new set of data from the Hipparcos survey (ESA 1997) allowed us to map the kinematics and variation of the stellar density close to the plane, resulting in a redetermination of the Oort limit (Crézé et al. 1998; Holmberg & Flynn 2004). A detailed bibliography of the most recent works about the K_{z} force, during the last decade, can be found in two recent publications (Garbari et al. 2012; Zhang et al. 2013), and in a review on the K_{z} problem and related questions by Read (2014).
This paper proceeds as follows. In Sect. 2, we present the sample selection, the determination of the vertical volume density of the RC stars, and their kinematics. Section 3 is devoted to the methods used to interpret the data. In Sect. 4, we present the results. Finally, conclusions are given and discussed in Sects. 5 and 6.
2. RC star selection, vertical density, and kinematics
2.1. The sample selection
A first sample is drawn from the fourth data release of the RAVE survey (Kordopatis et al. 2013) that contains stellar measurements of about 480 000 stars. Measurements for each star include radial velocity, effective temperature, gravity and metallicity. Crossidentification with the 2MASS and UCAC3 catalogues give us the corresponding photometry and proper motions. Our aim is to select RC stars towards the SGP within a 500 pcradius cylinder centred at the solar position, and with its main axis directed towards the SGP. The use of a cylinder for counts, instead of using a cone, allows us to increase the number of stars at low z, and to minimize the Malmquisttype bias. To minimize the consequence of interstellar extinction, we also restrict the selection to Galactic directions  b  > 22deg. Thus, taking into consideration the extinction A_{K} = 0.0116 and a scale height for the extinction of 90 pc (Groenewegen 2008), we neglect a systematic overestimation of distances that ranges between 0.6 and 1.5%. To include most of the RC stars and reject AGB stars, we apply as a colour selection J − K within the range from 0.5 to 0.8, based on 2MASS magnitudes (Fig. 1).
Fig. 1 Histograms of J − K colours for 12308 RAVE stars observed towards the South Galactic Pole (black continuous line). Red clump stars (red dashed), subgiants (blue dashdotted line), dwarfs (green dotted line) identified from their gravity. 
Fig. 2 Histogram of gravities for 12308 RAVE stars observed towards the South Galactic Pole. 
Based on highresolution spectroscopy of nearby Hipparcos stars, gravity measurements show that the red clump is defined by a restricted range of gravities from 1.8 to 2.8 (Soubiran et al. 2008; Valentini & Munari 2010). This is in agreement with expectations from stellar evolution models (Girardi et al. 2000), though they may suffer from uncertainties on the amount of mass loss during the first ascent of the giant branch. The histogram of gravities from the RAVE DR4 stars with J − K ∈ [ 0.5−0.8 ] shows a welldefined peak, allowing us to distinguish RC stars from subgiants and dwarfs (Fig. 2). We note a tail, containing a small fraction of stars with gravities smaller than 1.8, which might correspond to stars with lower absolute magnitudes along the red giant branch or along the AGB. Considering our colour selection interval, we do not expect such a significant number of these stars, seen neither among nearby Hipparcos stars nor in model predictions. Here, we believe that their presence mainly results from the uncertainties in the determination of RAVE gravities for giants. An indication for this is the small bias in this range of gravities shown in the Fig. 6 of Kordopatis et al. (2013), where RAVE and external gravities are compared. For this reason, we choose the gravity interval, 1.0 to 2.8, to define our RC sample. Using a sample by restricting log g within 1.8 and 2.8, reduces the sample by 20%. This does not modify our final conclusions concerning the measurement of the K_{z} force. It is a complementary indication that these stars, with RAVE log g ∈ [ 1.,1.8 ] are also RC stars.
Finally, we only select stars with a proper motion accuracy better than 4 mas y^{1} and better than 5 km s^{1} for the radial velocity, with signaltonoise ratio (S/N) greater than 20 in the RAVE spectra, and RV_{skycorrection}< 10. We reject stars without measurement of gravity or metallicity. In the case of multiple observations of the same star (about 10 percent of all measured stars), we keep the measurements with the highest S/N.
2.2. The RC stellar density
The RC stars are located in the HR diagram within restricted intervals of colour, absolute magnitude, and gravities. Using the last Hipparcos reduction (van Leeuwen 2007), Groenewegen (2008) determined the mean absolute magnitude of RC stars, and found M_{K} = −1.54 with a dispersion of 0.22 mag, but 0.15 just considering the stars with the most accurately measured parallaxes. As noted by Groenevegen, his results could be affected by the photometry of bright stars that are saturated in the 2MASS data. This was confirmed by Laney et al. (2012) who measured new IR photometry of 226 RC Hipparcos stars. Repeating the analysis, they obtained M_{K} = −1.613 ± 0.015. We note that the mean absolute Kmagnitude of RC stars does not depend on metallicity (Groenewegen 2008).
To define our main sample, we apply the supplementary selection criterion: ${\mathit{m}}_{\mathit{K}}\mathit{<}\mathrm{}\mathrm{1.613}\mathrm{+}\mathrm{5}\mathrm{log}\hspace{0.17em}\mathrm{\left(}\mathrm{500}\mathrm{\right)}\mathrm{}\mathrm{5}\mathrm{}\mathrm{5}\mathrm{log}\hspace{0.17em}\mathrm{\left(}\mathrm{cos}\mathrm{\right}\mathit{b}\mathrm{\left}\mathrm{\right)}\mathit{.}$(1)This criterion includes RC stars located within the 500 pcradius cylindric volume oriented towards the Galactic poles. The sample contains 9522 stars with J − K within [0.5,0.8], among which there are 5618 RC stars.
A second sample, extracted from the 2MASS catalogue, is used to determine the degree of completeness of our main RAVE sample and to quantify the selection function of RAVE observations. For this sample, the same criteria on J − K colour, apparent magnitudes (Eq. (1)), and Galactic directions are applied.
RAVE is a magnitudelimited survey of stars randomly selected in the southern celestial hemisphere. The original design was to only observe stars in the interval 9 <I< 12, but the actual selection function includes stars both brighter and fainter. For sufficiently small intervals of magnitudes and Galactic directions, RAVE stars are randomly selected^{1}. Hence, for these intervals, the fraction of RC stars in the sky can be directly estimated from the ratio of the number of RAVE RC stars to the number of observed RAVE stars.
We determine this ratio using RAVE counts in both north and south Galactic directions in order to improve the statistics at low z. The 2MASS counts are complete down to magnitude K = 14.3 (Skrutskie et al. 2006). Hence, we combine 2MASS counts, RAVE counts, and RAVE RC counts to estimate the total number of RC stars at any interval of K magnitudes and directions where RAVE observations exist.
In practice, since we cover a large range of Galactic latitudes and are interested in RC counts versus the vertical Galactic distance z, we do not determine the counts using the apparent magnitudes but using the corrected apparent magnitudes: ${\mathit{K}}_{\mathrm{C}}\mathrm{=}{\mathit{m}}_{\mathit{K}}\mathrm{+}\mathrm{5}\mathrm{log}\mathrm{\left(}\mathrm{sin}\mathrm{\right}\mathit{b}\mathrm{\left}\mathrm{\right)}\mathrm{=}{\mathit{M}}_{\mathit{K}}\mathrm{+}\mathrm{5}\mathrm{log}\mathrm{\left}\mathit{z}\mathrm{\right}\mathrm{}\mathrm{5}\mathit{,}$(2)which depends only on the absolute magnitude of the star and of its  z  position.
Fig. 3 2MASS star counts towards the South Galactic Pole versus the K_{c} magnitude (green circle) corrected for Galactic latitude (RAVE counts towards both poles are indicated by black crosses, and RAVE counts of red clump stars are denoted by red plus signs). 
Fig. 4 Vertical number density distribution of red clump stars towards the South Galactic Pole (black symbols: 200 pc binning, red symbols 50 pc binning). 
Figure 3 plots the 2MASS counts as a function of K_{c} magnitudes. Star counts from our main RAVE sample and from our RAVE subsample of RC stars are also plotted. From these three counts, we deduce the total number density of RC stars within the 2MASS sample, as a function of K_{c}, or equivalently as a function of the height  z  (Fig. 4). Error bars are determined from the three counts by using the statistical hypergeometric law.
Two sources of known bias are present but remain small in this analysis. The first known bias is the degree of homogeneity of the sample selections. Because of high S/N (K< 10), the accuracy of the various measured or used parameters remains high independent of the z distance. For instance, the median accuracy in J − K colours (within 0.5–0.8) is 0.03 from K = 6 to K = 10. Similarly, the mean S/N of the RAVE spectra used to determine the gravity remains high for RC stars at 2 kpc (K ~ 10): the mean S/N is 51 (rms 16). This implies that our selections and cuts remain homogeneous independent of the distance z.
A second effect is the Malmquist bias: this depends on σ_{M}, the dispersion of luminosity of the stellar candles, and on the variation of the density along the line of sight. In the case of a vertical exponential density law, ν ~ exp(−z/h), with h = 700 pc and σ_{M} = 0.2, at z = 1000 pc the bias on the estimated distances is +2% using a cone for the counts and is −0.7% using a cylinder. At z = 2000 pc the bias is +3% using a cone, and +1.2% using a cylinder. For the dynamical determination of the total mass perpendicular to the Galactic plane, we are interested in the density gradients, and so just in the variation of this bias: in this study, it is less than 1%. We note that with other tracers with an absolute magnitude dispersion of 0.5, the bias from star counts would be significantly larger: for cone counts, it is of the order of 5% at z = h and 11% at z = 3 h. This implies a systematic error of 6% on the resulting determination of the Galactic local surface mass density.
2.3. The RC star kinematics
We need to determine the vertical velocities of RC stars that combined with counts towards the Galactic poles will constrain the vertical potential at the solar position.
Radial velocities are obtained from our RAVE observations, proper motions from the UCAC3 catalogue and distances from the identification of RC stars (see Sects. 2.1 and 2.2). Radial velocities, proper motions, and distances of RAVE RC stars are converted in (u,v,w) velocities relative to the Sun, and in Galactic velocities, V_{R} − V_{⊙ ,R} and V_{z} − V_{⊙ ,z}, uncorrected for the solar motion, assuming R_{0} = 8.5 kpc.
The errors on the velocities are obtained from individual errors on proper motions and radial velocity, adopting a mean uncertainty on distances of 10% (Fig. 5). The median error on the V_{z} component is 2.4 km s^{1}.
Fig. 5 Distribution of errors in the vertical Galactic velocity for stars with  z  < 2000 pc (continuous black line), and for stars with 1300 <  z  < 2000 pc (dotted red line). 
Fig. 6 Vertical (black symbols) and radial (red symbols) velocity dispersions: σ_{Vz}, σ_{VR}. Mean vertical velocity V_{z} (black line). 
The mean vertical velocity is constant with z (Fig. 6). The velocity dispersions σ_{R} and σ_{z} are measured by applying a 3.5sigmaclipping to the V_{R}, V_{z} Galactic velocity components. The uncertainties on the dispersions are $\mathit{\sigma}\mathit{/}\sqrt{{\mathit{n}}_{\mathrm{\ast}}\mathrm{}\mathrm{1}}$. The vertical velocity dispersion σ_{Vz} rises up to 38 km s^{1} at 1 kpc and then remains nearly constant (Fig. 6).
The velocity ellipsoid tilt is null at z = 300 pc and reaches 8 ± 1deg at 1 kpc, pointing not far off the Galactic centre. This is in agreement with the finding by Siebert et al. (2008) and Pasetto et al. (2012a,b) based on a previous release of the RAVE survey. As discussed in Siebert et al. (2008) a bias on the measure of the tilt exists if no corrections are applied to consider the anisotropy in the errors of radial velocities and tangential velocities. This bias increases with distance (and with errors on tangential velocities), small at z = 1 kpc, it is about 7deg at z ~ 2 kpc. Unbiased measurements below 1 kc indicate that the tilt points not far from the Galactic centre. This is confirmed by recent measurements (Büdenbender et al. 2014) based on a large sample of SDSS/SEGUE G dwarfs for which the tilt increases with z and points in a direction close to the Galactic centre. This implies a correlation between radial and vertical motions, which we model using Stäckel potentials (see Sect. 3).
2.4. Metallicities
To improve the analysis of the vertical potential and the K_{z} force, we split our sample according to the metallicity in three subsamples delimited by the values [M/H] = −0.35 and −0.15. They contain 2182, 2558 and 2263 stars, respectively, of which 1440, 1741 and 1447 are RC stars. The K_{z} is mainly sensitive to the gradients of the density, and thus each subsample can probe more efficiently a different range of zdistances, the lowest metallicity sample probing the potential for the largest distances.
On the other side, the vertical potential determination is poorly determined at low z because of the lack of data below 300 pc. For this reason, we add a sample of Hipparcos and Elodie RC stars (300 stars) towards the North Galactic Pole. These stars were previously used to constrain the K_{z} force (Bienaymé et al. 2006; Soubiran et al. 2008). This small sample covers distances from 0 to 800 pc.
2.5. Secondary red clump
The tracer distance measurement is essential for the accuracy of the vertical potential determination.
The red clump giants span a large range in gravity depending on their metallicity: log g = 2.08 for the metalpoor, lowmass end and reaches up to log g = 3 for the highmass, metalrich red clump objects (Zhao et al. 2000).
Identification of clump stars from RAVE observations allow us to achieve a remarkable uncertainty of 7% in distances if we consider the dispersion of absolute K magnitude for Hipparcos clump stars with the best parallaxes (Groenewegen 2008), 10% otherwise. The existence of a secondary clump at slightly lower magnitude concerns higher mass stars (3 solar masses) for which the He burning does not operate in a degenerate gas (Girardi 1999). It may degrade the accuracy of our determination of distances. These stars are not very numerous compared to ordinary clump stars in the solar neighbourhood, they are younger and should not be present in older stellar populations with large velocity dispersions or z larger than about 300 pc.
This is in agreement with two recent findings. The first is the determination from asteroseismology (Stello et al. 2013) of the mass and of the main or secondary clump status for giants. The second is the comparison of histograms of asteroseismologic mass (Miglio et al. 2013) in two COROT fields at different low Galactic latitudes.
Finally, another source of possible bias is related to binarity, however, to modify significantly the apparent magnitude, binary systems consisting of giants are needed. Such systems are extremely rare, less than 1% (Nataf et al. 2012) since the life time on the red giant phase is extremely short, and the mass of each companion must not differ by more than 0.5% to have a binary consisting of two giants.
3. Methods and models
To measure the vertical potential at the solar Galactic position up to a vertical distance of 2 kpc, we use and adapt the centuryold method developed by Kapteyn (1922) and Oort (1932). Their method can be applied when the stellar oscillations through the Galactic plane remain smaller than ~1 kpc, so the vertical motions are approximately decoupled from the horizontal ones. Thus, below 1 kpc from the midplane, the problem becomes 1D dynamical model. In this case, the vertical distribution of stellar positions and vertical velocities f(z,w) can be written as the sum of isothermal components depending only on the vertical energy: $\begin{array}{ccc}\mathit{f}\mathrm{\left(}\mathit{z,w}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma}}_{\mathit{i}}\hspace{0.17em}\frac{{\mathit{c}}_{\mathit{i}}}{\sqrt{\mathrm{2}\mathit{\pi}}{\mathit{\sigma}}_{\mathit{i}}}\hspace{0.17em}\mathrm{exp}\left(\mathrm{}\frac{\mathrm{\Phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{}\mathrm{\Phi}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{w}}^{\mathrm{2}}}{{\mathit{\sigma}}_{\mathit{i}}^{\mathrm{2}}}\right)& & \end{array}$with the vertical density of the tracer stars as: $\begin{array}{ccc}\mathit{\nu}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{\Sigma}}_{\mathit{i}}\hspace{0.17em}{\mathit{c}}_{\mathit{i}}\hspace{0.17em}\mathrm{exp}\left(\mathrm{}\frac{\mathrm{\Phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{}\mathrm{\Phi}\mathrm{\left(}\mathrm{0}\mathrm{\right)}}{{\mathit{\sigma}}_{\mathit{i}}^{\mathrm{2}}}\right)\mathrm{\xb7}& & \end{array}$A general solution is provided in Appendix A. This solution may be written in a different, but equivalent form found by Garbari et al. (2012). Such a solution cannot be applied at higher z, where the coupling between radial and vertical motions must be considered in order to achieve an accurate measure of the vertical potential. For that purpose, Kuijken & Gilmore (1989) proposed analytic corrections that they applied to their sample of K dwarfs for the determination of the vertical force. Statler (1989) building exact stationary solutions with Stäckel potentials found corrections of the order of 10% at 1 kpc in comparison with a 1D model.
The closeness of the true Galaxy potential to a Stäckel potential is a very commonly used feature. For instance, Binney (2012) used this similarity to evaluate the three actions associated with any orbit, actions which he then uses as the arguments of his distribution functions. Here, we follow the approach of Statler (1989) and build a Galactic model with a Stäckel potential for which the HamiltonJacobi equation is fully separable, and a phasespace distribution function depending on the three straightforward integrals of motion in such a separable potential. By construction, the distribution function is stationary and solution of the collisionless Boltzman equation, and is used to model the vertical number density and vertical velocity dispersion of our stellar samples.
To model a realistic gravitational field within the solar neighbourhood, we follow the mass modelling proposed by Batsleer & Dejonghe (1994) and Famaey & Dejonghe (2003), using a combination of two Kuzmin & Kutuzov (1962) components, a disc and a halo (a detailed description of their properties can be found in Dejonghe & de Zeeuw 1988). We introduce four free parameters: the mass M and axis ratios ϵ of both components. These four free parameters allow us to constrain locally the scaleheight of the disc, its local surface density, and the local volume density of the extended component. Thus, the obtained modelled vertical distribution of the total volume density in the solar neighbourhood can represent any combination of a dark halo and a vertically extended stellar component. The fourth free parameter (the flattening of the halo component) is adjusted to impose a flat rotation curve over an extended range of Galactic radius. The K_{z} fit is thus just made on three parameters. We note that adjusting the flattening of the halo in this way does not modify its local density gradient, the density remaining nearly constant between z = 0 and 2 kpc, since the halo is never highly flattened. The halo flatness nevertheless affects the (poorly known) value of the corresponding circular velocity at the Sun’s radius.
Finally, with Stäckel potentials, we can very simply model the tilt of the velocity ellipsoid above the Galactic plane. The tilt orientation is fully determined by the positions ± z_{0} of the foci along the vertical axis, z_{0} defining a confocal ellipsoidal coordinate system. We set z_{0} = 2 kpc, in order to have a velocity ellipsoid at (R, z) = (8.5 kpc, 1 kpc) pointing close to the Galactic centre in agreement with observations (Siebert et al. 2008; Büdenbender et al. 2014).
3.1. KuzminKutuzov potentials
A detailed description of all characteristics of 3D Stäckel potentials can be found in de Zeeuw (1985). These potentials are easily tractable in confocal spheroidal coordinates. The prolate spheroidal coordinates (λ,θ,ν) are related to the cylindrical coordinates (r,θ,z) by: ${\mathit{r}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{(}\mathit{\lambda}\mathrm{+}\mathit{\alpha}\mathrm{)}\mathrm{(}\mathit{\nu}\mathrm{+}\mathit{\alpha}\mathrm{)}}{\mathit{\alpha}\mathrm{}\mathit{\gamma}}\mathrm{and}{\mathit{z}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{(}\mathit{\lambda}\mathrm{+}\mathit{\gamma}\mathrm{)}\mathrm{(}\mathit{\nu}\mathrm{+}\mathit{\gamma}\mathrm{)}}{\mathit{\gamma}\mathrm{}\mathit{\alpha}}\mathrm{\xb7}$(3)The shape of the coordinate surfaces is determined by α and γ while λ,ν satisfy −γ ≤ ν ≤ −α ≤ λ. Surfaces of constant λ are spheroids and those of constant ν are hyperboloids. They all share the same foci located on the z axis at ± z_{o} = ± (γ − α)^{1/2}.
A general Stäckel potential takes the form: $\mathrm{\Phi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\mathrm{}\frac{\mathit{h}\mathrm{\left(}\mathit{\lambda}\mathrm{\right)}\mathrm{}\mathit{h}\mathrm{\left(}\mathit{\nu}\mathrm{\right)}}{\mathit{\lambda}\mathrm{}\mathit{\nu}}$(4)where h is an arbitrary function.
Here, we define the class of KuzminKutuzov potentials (Dejonghe & de Zeeuw 1988), writing $\mathit{h}\mathrm{\left(}\mathit{\tau}\mathrm{\right)}\mathrm{=}\mathit{GM}\sqrt{\mathit{\tau}\mathrm{+}\mathit{q}}$ (with the condition q ≥ γ): ${\mathrm{\Phi}}_{\mathit{KK,q}}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\mathrm{}\frac{\mathit{GM}}{\sqrt{\mathit{\lambda}\mathrm{+}\mathit{q}}\mathrm{+}\sqrt{\mathit{\nu}\mathrm{+}\mathit{q}}}\mathrm{\xb7}$(5)The corresponding isodensity surfaces from Poisson equation are flattened oblate spheroids. Increasing ϵ^{2} = (q − α)/(q − γ) flattens the spheroids. The key is that adding multiple Stäckel potentials of this type still gives a Stäckel potential as long as the focal distance $\sqrt{\mathit{\gamma}\mathrm{}\mathit{\alpha}}$ remains the same.
Besides the energy and the angular momentum, a third independent isolating integral of the motion exists and can be written as: $\begin{array}{ccc}{\mathit{I}}_{\mathrm{3}}\mathrm{=}\mathrm{\Psi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}& & \mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\frac{{\mathit{z}}^{\mathrm{2}}}{\mathit{\gamma}\mathrm{}\mathit{\alpha}}\mathrm{(}\mathit{r\u0307}\mathrm{2}\mathrm{+}\mathrm{(}\mathit{r}\dot{\mathit{\theta}\mathrm{)}}{\mathrm{2}}^{\mathrm{)}}\\ & & \mathrm{}\frac{\mathrm{1}}{\mathrm{2}}\left(\frac{{\mathit{r}}^{\mathrm{2}}}{\mathit{\gamma}\mathrm{}\mathit{\alpha}}\mathrm{+}\mathrm{1}\right)\mathit{z\u0307}\mathrm{2}\mathrm{+}\frac{\mathit{rz}\mathit{r\u0307}\mathit{z\u0307}}{\mathit{\gamma}\mathrm{}\mathit{\alpha}}\end{array}$(6)with $\mathrm{\Psi}\mathrm{\left(}\mathit{\lambda ,\nu}\mathrm{\right)}\mathrm{=}\frac{\mathrm{(}\mathit{\nu}\mathrm{+}\mathit{\gamma}\mathrm{)}\mathit{h}\mathrm{\left(}\mathit{\lambda}\mathrm{\right)}\mathrm{}\mathrm{(}\mathit{\lambda}\mathrm{+}\mathit{\gamma}\mathrm{)}\mathit{h}\mathrm{\left(}\mathit{\nu}\mathrm{\right)}}{\mathrm{(}\mathit{\gamma}\mathrm{}\mathit{\alpha}\mathrm{)}\mathrm{(}\mathit{\lambda}\mathrm{}\mathit{\nu}\mathrm{)}}\mathrm{\xb7}$(7)
3.2. The distribution function
To model the density and kinematics of our samples, we define a stationary distribution function that depends on three integrals of the motion. We use the 3D stellar disc distribution function of Bienaymé (1999), which has nearly a Schwarzschild distribution behaviour in the limit of small velocity dispersions. This distribution is a generalization of the Shu distribution function (Shu 1969) and also has a density that is nearly radially exponential.
The distribution function is
$\begin{array}{ccc}& & \mathit{f}\mathrm{\left(}\mathit{E,}{\mathit{L}}_{\mathit{z}}\mathit{,}{\mathit{I}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}\mathrm{\Omega}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}}{\mathrm{2}\mathit{\pi \kappa}\mathrm{\left(}{\mathit{R}}_{\mathrm{c}}\mathrm{\right)}}\frac{\mathrm{\Sigma}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}}{{\mathit{\sigma}}_{\mathit{r}}^{\mathrm{2}}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}}\mathrm{exp}\left[\mathrm{}\frac{\mathit{E}\mathrm{}{\mathit{E}}_{\mathrm{circ}}}{{\mathit{\sigma}}_{\mathit{r}}^{\mathrm{2}}}\right]\\ & & \mathrm{\times}\hspace{0.17em}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi}}}\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{z}}\mathrm{\left(}{\mathit{L}}_{\mathit{z}}\mathrm{\right)}}\mathrm{exp}\begin{array}{c}\mathrm{\u23a7}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23a8}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23a9}\end{array}\mathrm{}{\left(\frac{{\mathit{R}}_{\mathrm{c}}\mathrm{(}{\mathit{L}}_{\mathit{z}}{\mathrm{)}}^{\mathrm{2}}}{{\mathit{z}}_{\mathit{o}}^{\mathrm{2}}}\mathrm{+}\mathrm{1}\right)}^{1}\left(\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}\mathrm{}\frac{\mathrm{1}}{{\mathit{\sigma}}_{\mathit{r}}^{\mathrm{2}}}\right){\mathit{I}}_{\mathrm{3}}\begin{array}{c}\mathrm{\u23ab}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ac}\\ \mathrm{\u23aa}\\ \mathrm{\u23aa}\\ \mathrm{\u23ad}\end{array}\end{array}$(8)with R_{c}(L_{z}) the radius of the circular orbit that has the angular momentum L_{z}, Ω is the angular velocity, κ is the epicyclic frequency, and E_{circ} is the energy of a circular orbiting star at radius R_{c}.
For sufficiently small velocity dispersions, the number density distribution, Σ(L_{z}) = Σ_{0}exp(−R_{c}/R_{ν}), is close to Σ(R) = Σ_{0}exp(−R/R_{ν}).
We set σ_{r,z}(L_{z}) = σ_{0,r,z}exp(−R_{c}/R_{σr,z}) and the velocity dispersions are close to σ_{r,z}(R) = σ_{0;r,z}exp(−R/R_{σ}).
Local number density and dispersion, Σ_{0}, σ_{0}, are constants, R_{ν} is (close to) the scale length of the number density distribution, and R_{σ} is (close to) the scale length of the velocity dispersions.
This distribution function is very similar to that proposed by Statler (1989), but here, generalized to the cases where R ≠ R_{0}. We also reduced the number of free parameters to have a form closer to the Shu (1969) distribution function. The corresponding velocity distribution is not far from a 3D Gaussian. When z = 0, the corresponding density varies nearly exponentially in an extended domain of a few kpc around the Sun. We note that if z ≠ 0, the density above the plane also depends on the vertical potential; thus the density may vary exponentially at any z with a supplementary condition, as for instance R_{σ} = 2R_{ρ} (see also Eqs. (6), (7) of van der Kruit & Freeman 2011). The velocity dispersions also decrease exponentially. This distribution function had been previously used for a dynamically consistent analysis of the kinematics of Hipparcos stars (Bienaymé 1999).
For the modelling of the distribution function of our RC star samples, we fix the scale lengths for the radial density and for the velocity dispersions: R_{ν} = 2.5 kpc and R_{σ} = 9 kpc. This is to be compared with R_{ν} ~ 2.2 and 2.8 kpc for the thin and thick discs (CabreraLavers et al. 2005; Jurić et al. 2008; Chang et al. 2011; Polido et al. 2013; Robin et al., in prep.) and R_{σ2} = 4.4 or 5.6 kpc, (R_{σ} = 8.8 or 11.1 kpc), respectively by Lewis & Freeman (1989) and Ojha et al. (1996). We also set σ_{r}/σ_{z} = 2 for the thin disc stars (σ_{z}< 27 km s^{1}) and 1.5 otherwise, in agreement with the observed properties of our RC sample.
3.3. The corrected bias
The introduction of a 3D model allows us to correct different effects relative to a 1D vertical model. The first effect is the velocity ellipsoid tilt at large z that increases the observed vertical velocity dispersion. The second effect is related to the vertical bending of the stellar orbits. For stars seen at z = 1 or 2 kpc, their mean Galactic radius <R>, when they cross the Galactic plane, is larger than R_{0}, a position where the stellar density is lower than at the radius R_{0}. These two effects lead to an overestimate of K_{z} in a 1D model. A third effect is also related to the bending of orbits and to the radial gradient of σ_{z}, lowering σ_{z} towards the pole. This effect leads to an underestimate of the K_{z} force in a 1D model.
Here, the modelling with a locally valid Stäckel potential allows us to correct in a dynamically consistent way the bias resulting from the coupling of vertical and horizontal stellar motions. However, this is obtained at the expense of supplementary parameters: the tilt orientation of the velocity ellipsoid, the radial gradients of the stellar density, and of the kinematics.
3.4. The adjustment procedure
We determine the parameters of the vertical potential by fitting the observed moments, density ν(z) and vertical velocity dispersions σ_{zz}(z) for each of the three metallicity samples; moments are computed from the distribution function Eq. (8). We minimize the difference between observed and modelled quantities using the χ^{2}: $\begin{array}{ccc}{\mathit{\chi}}^{\mathrm{2}}\mathrm{=}{\mathit{\chi}}_{\mathit{\nu}}^{\mathrm{2}}\mathrm{+}{\mathit{\chi}}_{\mathit{\sigma}}^{\mathrm{2}}\mathrm{+}\mathit{q}{\left(\hspace{0.17em}{\mathit{V}}_{\mathrm{c}}\mathrm{(}\mathrm{9.5}\mathrm{kpc}\mathrm{)}\mathrm{}{\mathit{V}}_{\mathrm{c}}\mathrm{(}\mathrm{7.5}\mathrm{kpc}\mathrm{)}\hspace{0.17em}\right)}^{\mathrm{2}}& & \end{array}$with $\begin{array}{ccc}{\mathit{\chi}}_{\mathit{\nu}}^{\mathrm{2}}& \mathrm{=}& {\mathrm{\Sigma}}_{\mathit{i}}\frac{{\left({\mathit{\nu}}_{\mathrm{mod}\mathit{,i}}\mathrm{}{\mathit{\nu}}_{\mathrm{obs}\mathit{,i}}\right)}^{\mathrm{2}}}{{\mathit{\u03f5}}_{\mathit{\nu ,i}}^{\mathrm{2}}}\\ {\mathit{\chi}}_{\mathit{\sigma}}^{\mathrm{2}}& \mathrm{=}& {\mathrm{\Sigma}}_{\mathit{i}}\frac{{\left({\mathit{\sigma}}_{\mathrm{mod}\mathit{,i}}\mathrm{}{\mathit{\sigma}}_{\mathrm{obs}\mathit{,i}}\right)}^{\mathrm{2}}}{{\mathit{\u03f5}}_{\mathit{\sigma ,i}}^{\mathrm{2}}}\hspace{0.17em}\mathit{,}\end{array}$where the ϵ_{i} are the corresponding uncertainties.
For the density, the bin size is 100 pc between 300 pc and ~2000 pc. We do not consider our density estimates below 300 pc where the completeness is difficult to determine. For the velocity dispersions, the bin size is 100 pc between 200 pc and 1200 pc and 200 pc beyond.
The last r.h.s. term within the χ^{2} expression is introduced to impose a potential with a nearly flat rotation curve. The q factor does not need to be large and the contribution of the corresponding term to the χ^{2} is small, because the slope of the rotation curve is uncorrelated to the other parameters. The adjusted rotation curve is flat and varies by 1 km s^{1} between R = 7 and 15 kpc.
We use two quasiisothermal components (Eq. (8)) to model the stellar sample distribution function, the density ν(z) and dispersion σ_{z}(z) for each metallicity sample.
We determine the χ^{2} minimum for the parameters of the potential using the MINUIT software (James 2004) that allows us to look for possible multiple minima and to obtain a first estimate of the variancecovariance matrix. To compute the posterior probability distribution function (PDF), we consider the likelihood ℒ = exp(−χ^{2}/ 2) and apply a Markov chain MonteCarlo using the MetropolisHastings Algorithm (ForemanMackey et al. 2013). From this, we determine the marginal PDFs.
Fig. 7 Vertical number density (top) and velocity dispersions (bottom), data (symbols with error bars) and model (lines) for the three RAVE metallicity samples and the NGPHipparcos sample: low metallicity [ M/H ] ≤ −0.35 (black), intermediate −0.35 < [ M/H ] ≤ −0.15 (red), high −0.15 < [ M/H ] (green), NGP (blue). Dotted lines are the extrapolated model not fitted to data. 
Fig. 8 Vertical number density (top) and velocity dispersions (bottom) for the full RAVE RC sample (symbols with error bars) and model (continuous line). 
Figure 7 shows the binned data and the best fit model to the density and dispersion versus z for the three metallicity RAVE samples and for the HipparcosNGP sample. For information, we also show this bestfit model and the binned data for density and vertical velocity dispersions without metallicity splitting (Fig. 8).
4. Results
4.1. The K_{z} force
In the previous sections, we have described the stellar samples used to probe the vertical force towards the Galactic pole. Our method does not fundamentally differ from the centuryold pioneering work of Oort (1932). The modelled quantity that is fitted to the observations is the vertical potential within the interval of distances probed by our stellar tracers. The vertical force and the total mass density distributions are deduced from the first and second zderivatives of this potential. It is known that the K_{z} is an illconditioned problem and that without a filtering of data or smoothing assumptions about the shape of the potential, the derivatives would be dominated by the noise and fluctuations of data. Here, the actual smoothing assumption for the local potential is given by two KuzminKutuzov components (Batsleer & Dejonghe 1994; Famaey & Dejonghe 2003) to mimic the potential of a disc and a spheroid (and also a flat rotation curve). This is partly equivalent to the threeparameter modelling introduced by Kuijken & Gilmore (1989) to describe the total Galactic disc mass surface density. Our model goes beyond the planeparallel approximation and allows us to describe the distribution function f(z,w) beyond 1 kpc up to 2 kpc. The vertical potential is modelled by a disc, its local density and thickness, and by the local density of an extended component that represents the dark matter halo.
Fig. 9 K_{z} force and 1σ error intervals. 
The resulting K_{z} at R_{0} from the fit of RAVE data (Fig. 7) is shown in Fig. 9. We give the results in terms of K_{z}/ 2πG for visibility and comparison with other studies, even though this should not be confused with the true surface density.
The obtained K_{z} is less constrained at low z where there are no RAVE data, but only the smaller sample of Hipparcos and Elodie stars. The K_{z} force at 350 pc is found to be 44.2$\begin{array}{c}\mathrm{+}\mathrm{2.3}\\ 2.9\end{array}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot}}$ pc^{2} in agreement with the Korchagin et al. (2003) determination 42 ± 6 M_{⊙} pc^{2}, based on a sample of 1500 K giants from Hipparcos observations, mainly first ascent giants rather than clump giants.
We also note that our K_{z} determination near z = 0 allows us to deduce the Oort limit ρ_{dyn}(z = 0) = 0.0911 ± 0.0059 M_{⊙} pc^{3}, intermediate value between the Crézé et al. (1998; 0.076 ± 0.015 M_{⊙} pc^{3}) and Holmberg & Flynn (2004; 0.102 ± 0.010 M_{⊙} pc^{3}) determinations. These two studies are based on the same Hipparcos data, below z = 125 pc, but with different assumptions regarding the shape of the potential. Most likely, this explains the difference between both results, which differ by just a 1σ error.
Our K_{z} force determination from 0 to 1 kpc, is similar to recent studies, but in our case, with more free parameters and without limiting assumptions on the baryonic local surface density or on the dark matter local volume density.
Since the quasitotality of the ordinary matter resides below z = 1 kpc, the mass density beyond 1 kpc is dominated by the dark matter. This results that our K_{z} measurement gives direct access to the dark matter density between 1 and 2 kpc above the Galactic plane.
The K_{z} force at intermediate and high z distances are: $\begin{array}{ccc}{\mathit{K}}_{\mathit{z}}\mathrm{(}\mathrm{1}\mathrm{kpc}\mathrm{)}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{\pi}\mathrm{G}\mathrm{\right)}& \mathrm{=}& \mathrm{68.5}\mathrm{\pm}\mathrm{1.0}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{2}\\ & \mathrm{=}& \mathrm{1852}\mathrm{\pm}\mathrm{27}{\mathrm{km}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{s}}^{2}\hspace{0.17em}{\mathrm{kpc}}^{1}\mathit{,}\end{array}$and $\begin{array}{ccc}{\mathit{K}}_{\mathit{z}}\mathrm{(}\mathrm{2}\mathrm{kpc}\mathrm{)}\mathit{/}\mathrm{\left(}\mathrm{2}\mathit{\pi}\mathrm{G}\mathrm{\right)}& \mathrm{=}& \mathrm{96.9}\mathrm{\pm}\mathrm{2.2}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{2}\\ & \mathrm{=}& \mathrm{2619}\mathrm{\pm}\mathrm{59}{\mathrm{k}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{s}}^{2}\hspace{0.17em}\mathrm{kpc}}^{1}\mathit{.}\end{array}$
4.2. Dependence on fixed model parameters
The vertical tilt of the velocity ellipsoid is fixed in our models to point close to the Galactic centre through the choice of our foci (tilt of ~13 deg at z = 2 kpc, in accordance with Büdenbender et al. 2014). Unbiased measurements indicate that the tilt indeed points close to the Galactic centre, and our method is not affected by the existence of a bias in the tilt that one would actually measure with RAVE data at large heights (Sect. 2.3) because we directly fit individual velocities (Sect. 3.4) and not the global shape of the velocity ellipsoid.
The two other fixed parameters, nonexistent in traditional analysis assuming separation of vertical and radial motions, are the radial scale lengths of the number density of tracer stars and of the velocity dispersions.
The radial scale length of the stellar sample, R_{ν}, modifies by less than 1% the K_{z} between z= 0 to 2 kpc, when R_{ν} is varying from 1.8 kpc to 3.5 kpc. Increasing the kinematics scale length does not modify the K_{z} determination. Only decreasing the scale length R_{σ} to 7 or 5 kpc (or R_{σ2} to 3.5 or 2.5 kpc) increases the K_{z} force at 2 kpc by 5% and 11%, respectively. This implies that our determination of the surface mass density of disc would be be increased by 5% or 10%, and the local DM density by 5% or 14%. However, these small kinematic scale lengths are excluded by existing observations.
4.3. The vertical mass density
To be able to estimate the vertical mass density distribution of Galactic components, the K_{z} determination is not completely sufficient, and we must also know the 3D shape of the baryonic and dark matter components. Here, the fourparameter Stäckel potential we fitted should be considered simply as a way to estimate the K_{z} force itself, but the relative contribution of the baryonic disc and halo to this force can be more reliably dealt with a posteriori by representing the baryonic mass component with a double exponential law ρ(R,z) ~ exp(−R/H_{ρ})exp(− z  /h_{z}), whose mass is assumed to be proportional to the stellar discs, the dominant baryonic mass component. Recent analyses of Galactic star counts with accurate and detailed modelling of the luminosity functions converge towards a short scale length, 2.1 to 2.3 kpc for the thin disc (see for instance Robin et al., in prep.) and between 2.8 to 3.2 for the thick disc, while the vertical density of the stellar disc population is very close to an exponential (Fig. 10).
Fig. 10 Vertical mass density of the stellar thin discs from the Besançon model (red line) fitted with an exponential of scale height h = 229 pc (black line). 
For a given mass surface density of the disc at R_{0}, both parameters H_{ρ} and h_{z} can modify the vertical force. Decreasing the scale length H_{ρ} increases the vertical force of the disc. In this case, to fit the observed K_{z}, the surface mass density of the disc must be decreased (and to fit the K_{z} at larger z, the dark matter mass density ρ_{DM} must be increased).
Here, we will consider as reference values H_{ρ} = 2200 pc, h_{z} = 300 pc and R_{0} = 8500 pc for the double exponential disc.
We assume that the dark halo component is spherical. Its radial mass density is defined to exactly complement the doubleexponential disc component in order that our model of the Galactic rotation curve is strictly flat.
Two free quantities remain: the mass or the local density of each component. We adjust the local density of the baryonic and of the dark matter to fit in a least square sense the observed K_{z}. The resulting total vertical mass density distribution is shown in Fig. 11 with the DM halo and disc decomposition. The plotted surface densities are the integrated mass volume density between −z and z. We obtain for the total surface density of the disc component at the solar position: $\begin{array}{ccc}{\mathrm{\Sigma}}_{\mathrm{disc}}\mathrm{\left(}{\mathit{R}}_{\mathrm{0}}\mathrm{\right)}\mathrm{=}\mathrm{44.4}\mathrm{\pm}\mathrm{4.1}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{2}\mathit{.}& & \end{array}$The probability distribution function of total local volume mass density of the dark matter component is plotted in Fig. 12. This yields $\begin{array}{ccc}{\mathit{\rho}}_{\mathrm{DM}}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{0.0143}\mathrm{\pm}\mathrm{0.0011}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{3}\mathrm{=}\mathrm{0.54}\mathrm{\pm}\mathrm{0.004}\hspace{0.17em}\mathrm{Gev}\hspace{0.17em}\mathrm{c}{\mathrm{m}}^{3}\mathit{.}& & \end{array}$In the case of the baryonic and the total (baryonic+DM) local volume densities, the Oort limit, we obtain
$\begin{array}{ccc}& & {\mathit{\rho}}_{\mathrm{baryons}}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{0.077}\mathrm{\pm}\mathrm{0.007}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{3}\mathit{,}\\ & & {\mathit{\rho}}_{\mathrm{total}}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathrm{)}\mathrm{=}\mathrm{0.091}\mathrm{\pm}\mathrm{0.0056}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{3}\mathit{,}\end{array}$and Fig. 13 plots their respective probability distribution function.
Fig. 11 Total surface mass density (black continous line) $\mathrm{\Sigma}\mathrm{(}\mathit{<}\mathrm{}\mathit{z}\mathrm{\left}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int}}}_{\mathrm{}\mathit{z}}^{\mathit{z}}{\mathit{\rho}}_{\mathrm{tot}}\mathrm{d}\mathit{z}$ at R_{0} split in a DM spherical component (red dashed line) and in a baryonic double exponential disc (blue dashdotted line). Dotted lines are 1σ error intervals. 
Fig. 12 Probability distribution function for the local dark matter density ρ_{DM}(z = 0). 
Fig. 13 Probability distribution functions for the total local mass density, i.e. the Oort limit ρ_{tot}(z = 0) (continous line), and the baryonic mass density ρ_{bar}(z = 0) (dotted line). 
Table 1 shows the results for ρ_{DM} and Σ_{D} in the case of some other disc parameter values (results are based on a χ^{2} minimization and are slightly different from these resulting from MCMCs). The change of the bestfit solution with the chosen method could result from the imperfect adequacy of the model. This would indicate a bias; this dependency of the results on the adopted methods has been also shown in a different context by Polido et al. (2013) analysing Galactic star counts.
Diskhalo parameters reproducing the observed K_{z} vertical force.
We can note that, remarkably, the dark matter density is nearly insensitive to the model parameters H_{ρ}, h_{z}, R_{0}. This results from the fact that the dark matter model has a nearly constant density in the range z = 0 to 2 kpc and depends quite exclusively on the difference on the K_{z} force between 1 and 2 kpc.
5. Discussion
The last RAVE data release (DR4, Kordopatis et al. 2013) allows us to probe the vertical density distribution of RC stars to a distance of 2 kpc from the Galactic plane, and also to determine their vertical kinematics and metallicity. This provides a highly accurate sample for the study of the vertical force perpendicular to the Galactic plane. About 5000 RC stars are used, permitting us to relax some of the traditional assumptions. Specifically, because of the large range of vertical distances probed up to 2 kpc, we can separate the contribution to the vertical force due to the halo and the disc.
The K_{z} problem is, in principle, illconditioned. The potential is the fitted quantity, and to determine the total mass vertical density, we evaluate the second derivative of the potential. To avoid arbitrary fluctuations resulting from the finite size of our stellar tracers, it is necessary to smooth the potential to recover a realistic vertical density. Here, this is done by assuming the local vertical potential shape as a Stäckel combination of a disc and a spheroidal halo. This is accomplished with a threeparameter model constraining the local densities of a disc and of a spheroidal halo, and the thickness of the disc. A fourth parameter imposes a flat rotation curve. The flatness of the rotation curve does not directly impact the K_{z} determination, but allows us a more realistic dynamical selfconsistent distribution function modelling the density and velocity distribution of tracer stars.
5.1. The K_{z} force measurement
At z = 1 kpc, our K_{z} measurement is similar and in agreement with the two last decades determinations of the K_{z} force based on data from stars up to z ~ 1 kpc: K_{z}(1 kpc/2π G) ~ 68 M_{⊙} pc^{2} (for recent determinations, see Garbari et al. 2012; Bovy & Rix 2013; Zhang et al. 2013). Between z = 1 and 2 kpc, we find that the scaleheight of the dark halo is significantly larger than the domain probed with our sample, and that the density is nearly constant below z = 2 kpc.
The local density of the halo is ρ_{DM}(z = 0) = 0.0143 ± 0.0011 M_{⊙} pc^{3} = 0.542 ± 0.042 Gev cm^{3}. We also determine the total surface density of the disc components Σ_{bar} = 44.4 ± 4.1 M_{⊙} pc^{2} that is in agreement with the Flynn et al. (2006) determination. The halo volume density and disc surface density are determined nearly independently, with the consequence that their resulting uncertainties are small. Nevertheless, systematic uncertainties, e.g. due to the choice of smoothing of the potential, are, in fact, larger.
Most of the recent determinations claimed estimates of the local halo density to be of the order of ρ_{DM}(z = 0) = 0.006−0.008 M_{⊙} pc^{3} (see the compilation in Table 4 by Read 2014) However, it must be noted that in these previous K_{z} measurements, the z extension of the stellar tracers was limited to distances of 0.8 to 1.1 kpc. For that reason, it was not possible to accurately separate the respective contributions from the visible discs (stellar and ISM discs) and those contributions from the dark matter halo. Then independent estimates of the visible matter from star counts and kinematics were adopted. As mentioned in the review by Read (2014), these determinations of the local dark matter density were based on an assumed total surface disc density for the visible matter. In one study, Kuijken & Gilmore (1989) assumed the value of the density of the dark matter and deduced the surface mass density of the disc.
5.2. The local volume mass density determination
To move from the vertical force to the vertical density distribution and to the decomposition in contributions from disc and from halo components implies we know nonlocal Galactic characteristics as the scale length of the disc, the Galactic solar radius R_{0}, or the thickness of the visible matter. Unfortunately, many fundamental characteristics of our Galaxy’s structure remain inaccurately measured, such as the distance of the Sun to the Galactic centre, or are even subject to contradictory determinations, such as the amplitude and shape of the Galactic rotation curve. For this reason, we list the consequences of our findings of high local density of dark matter, according to different hypotheses and to results recently published concerning the Galactic gravitational potential.
First, the large scale and 3D properties of the dark matter distribution are mainly established from the knowledge of the Galactic rotation velocity curve or from orbits of streams through the Galactic halo. Other constraints exist locally, for instance the Galactic escape velocity in the solar neighbourhood (Smith et al. 2007; Piffl et al. 2014a) or the determination of Oort’s constants (Olling & Merrifield 1998).
The flatness of the Galactic rotation velocity curve at R_{0} is consistent with the observations of external disc galaxies of similar type. The flatness is supported by the two recent determinations of the Galactic circular velocity curve from parallaxes (Reid et al. 2014) or from spectrophotometric distances (Bovy et al. 2012) that favour a flat rotation curve in an extended interval of radii around R_{0}.
Now, considering such a flat rotation curve with V_{c} ~ 220 km s^{1}, similar to the recent determination from the APOGEE project (Bovy et al. 2012), and a spherical dark matter halo (frequently called the “standard dark Galactic halo model” within the astroparticles literature), Salucci et al. (2010) and others built Galactic models and estimated ρ_{DM} ~ 0.006−0.01 M_{⊙} pc^{3} in the solar neighbourhood (see Table 4 from Read 2014).
If we consider the Galactic rotation curve (Reid et al. 2014) deduced from data of the BeSSeL project (Brunthaler et al. 2011), which leads to a higher circular velocity curve of about V_{c} ~ 240 km s^{1} (McMillan 2011; Bobylev & Bajkova 2013; Reid et al. 2014), and also assuming a spherical dark matter halo, McMillan (2011) obtained ρ_{DM} = 0.40 ± 0.04 Gev cm^{3}(=0.010 M_{⊙}/pc^{3}) a relatively low value still.
Concerning the local slope of the rotation curve, some of the strongest evidence comes from the determination of Oort’s constants by Olling & Dehnen (2003), based on Hipparcos data (see also Mignard 2000). Their determination is based on the oldest red giant stars, over a large domain of 2 kpc radius around the Sun. Their careful analysis considers the necessary corrections for the bias due to the extinction and also to the mode mixing from the solar motion. They found out that the rotation curve is flat over ±2 kpc from the Sun’s Galactic radius, varying by less than (A + B)/2 = dV_{c}/ dR = +0.5 ± 0.8 km s^{1} kpc^{1} (and dlnV_{c}/ dlnR = +0.02). Their value of the slope of velocity curve introduces a small additional nonzero contribution to the local mass density from the Poisson equation: $\begin{array}{ccc}\mathrm{(}\mathrm{4}\mathit{\pi GR}{\mathrm{)}}^{1}\frac{\mathit{\partial}{\mathit{V}}_{\mathrm{c}}^{\mathrm{2}}}{\mathit{\partial R}}\mathrm{=}\frac{{\mathit{B}}^{\mathrm{2}}\mathrm{}{\mathit{A}}^{\mathrm{2}}}{\mathrm{2}\mathit{\pi G}}\mathrm{=}\mathrm{+}\mathrm{0.0012}\mathrm{\pm}\mathrm{0.0019}{\mathit{M}}_{\mathrm{\odot}}\hspace{0.17em}{\mathrm{pc}}^{3}\mathit{.}& & \end{array}$For high values of the local dark matter density such as those we obtained here, we can mention the work by Burch & Cowsik (2013) who built a global dynamical Galactic model to constrain ρ_{DM} at the Galactic centre and at R_{0} from published observations of the circular velocity curve. They plot the K_{z} force that can be directly compared with our measurement. One of their K_{z} models (Fig. 6b) is in agreement with our determination within z = 0 to 2 kpc and they obtained ρ_{DM}(z = 0) = 0.015 M_{⊙} pc^{3}. Unfortunately, in their model they adopt a rapidly rising rotation curve at R_{0} that we judge very unlikely. This rising rotation curve, with dlnV_{c}/dlnR ~ 0.16 explains their high value obtained for ρ_{DM}. In their model, if the rotation curve was flat, the local density would probably be about 0.006 M_{⊙}pc^{3}.
In a recent study, using tracer stars with z between 1 to 2 kpc, Smith et al. (2012) estimated an order of magnitude for ρ_{DM}(z = 0) = 0.015 M_{⊙} pc^{3}, similar to our finding. However, because of their crude modelling of the potential and the lack of a clear definition of the selection function, they decided not to give error bars for their estimate. It remains remarkable that this unique previous analysis of distant tracer stars agrees with our finding. Finally, Piffl et al. (2014b) also recently determined a similar value as ours with the RAVE data, for a halo flattening of 0.8.
5.3. Global galactic properties
When we consider our model with ρ_{DM} = 0.014 M_{⊙} pc^{3}, a flat rotation curve and a spherical halo, this implies V_{c} = 267 km s^{1}, which is a too large value compared to the recent determinations of the Galactic rotation curve. The simplest way to reconcile our local determination of the dark matter density with the admitted flat rotation curves based on direct observations consists in flattening the dark matter halo. An axis ratio of the order of 0.67 is necessary, in the case of V_{c} = 220 km s^{1} as determined by Bovy et al. (2012). This significant flattening is not in agreement with the results of Galactic cosmological numerical simulations (Macciò et al. 2007) for which a mean flattening of the order of 0.8 is expected. Moreover, Pillepich et al. (2014) analysed simulations including dissipational gas physics and obtained much rounder halo with q = 0.99, instead of q = 0.53 in the case of DMonly numerical simulations.
In fact, combining the circular velocity from the BeSSeL project V_{c} ~ 240 km s^{1} (McMillan 2011; Bobylev & Bajkova 2013; Reid et al. 2014), with a nearly flat rotation curve at the Sun position, and a flattening of 0.8, leads to our estimated value of the local dark matter density, in accordance with Piffl et al. (2014b). This value of the circular velocity makes the Milky Way a clear outlier from the TullyFisher relation (Holmberg et al. 2006; Hammer et al. 2012).
Another plausible explanation of a high local DM density comes from the cosmological numerical simulations by Read et al. (2009, Fig. 9) and Pillepich et al. (2014; see also Ling 2010, for DM detection implications). They showed that in the case of a disc galaxy already in place at high redshift, the later accretion of galaxy satellites create a slowly rotating very thick disc or flattened spheroidal component of dark matter. This dark component results from the accretion of the dark component of each accreted satellite. Due to the history of accretion, this accreted DM component has a high angular momentum, with kinematical properties intermediate between the stellar disc and a nonrotating spherical dark halo. Its detailed structure depends on the details of the accretion history, and the halo mass depends on the unknown number and mass of accreted satellites. The local contribution of this accreted DM component could be 25 to 150 percent the density of the primordial and nearly spherical DM component. In the case of a small scale height of the order of 2–3 kpc, its vertical density is quickly decreasing. In such a case, our modelling of the K_{z} force probably does not include a sufficient number of free parameters to accurately model the shape of the K_{z} force between 1 and 2 kpc. We might suspect that our modelling forces a nearly linear rise of the K_{z} between 1 and 2 kpc. A supplementary (thick disc) KuzminKutuzov component could be added (Famaey & Dejonghe 2003) to model more precisely the vertical dark matter density and potential, but it is not clear that the size of our sample will allow us to discriminate between models with so many parameters.
A similar disc of “phantom” dark matter (from the point of view of a Newtonian observer) is predicted (Bienaymé et al. 2009) by the MOND effective theory. This can be a very similar effect to that observed in numerical simulations of accretion of satellite galaxies. For a baryonic model like that assumed here, a vertical force K_{z}/(2πG) ~ 90 M_{⊙} pc^{2} is predicted at z = 2 kpc. Nevertheless, a somewhat high value ~75 M_{⊙} pc^{2} is rather predicted at z = 1 kpc. For the same reason mentioned previously, our two KuzminKotozov components modelling of the vertical force is not adequate to correctly reproduce the K_{z} force in the case of Mondian model that shows a significant bending of the vertical force law between z = 1 and 2 kpc.
5.4. A Galactic DM halo with core?
Now, if we consider the consistency of our local K_{z} determination with ΛCDM, we first note that our result is nearly independent from any assumptions on the form of the Galactic potential, for instance in the central regions of the Galaxy. The decomposition of the K_{z} force in dark and visible contributions requires additional information, for instance, the disc scale length of the visible matter. Previous determinations indicated that there was less dark matter mass than predicted by cosmological simulations within R_{0} (Navarro & Steinmetz 2000; Binney et al. 2000; Famaey & Binney 2005; Abadi et al. 2010), whilst our determination would imply that there is more dark matter mass inside R_{0}. This has implications for the mass concentration and can be tested relative to the massconcentration relation reported by Macciò et al. (2008). This point is dicussed in detail in Piffl et al. (2014b) and they concluded that modifying the Galactic halo profile by taking NFW adiabatic contraction into account, they can obtain an agreement with simulations in a ΛCDM universe. We also remark that the recent Galactic DM halo modelling by Nesti & Salucci (2013) constrained with inner terminal velocities, MASER observations, and stellar halo velocity dispersions gives a high local DM mass density consistent with our finding. Their modelling favoured a cored profile (R_{H} ~ 10 kpc) of the DM halo (see also Bissantz et al. 2003).
6. Conclusion
We have established that a significant amount of dark matter resides close the Galactic disc: the local dark matter mass density is ρ_{DM}(z = 0) = 0.0143 ± 0.0011 M_{⊙} pc^{3}. We have independently determined the local DM density and the baryonic disc surface density at the solar Galactic radius R_{0}. The large size of our sample leads to small statistical errors, but it is clear that systematic errors could also arise from neglected elements in our modelling. For instance, one aspect of dynamical modelling can be questioned: resonant orbits (from vertical relative to horizontal motions) are numerous at z beyond 1 kpc and are not modelled by Stäckel potentials, or with a simple torus fitting. Also, the nonaxisymmetric effects due to spiral arms (Faure et al. 2014; Debattista 2014) or nonequilibrium features generated by the potential interaction of satellites with the Galactic disc (Gómez et al. 2013) could also bias the result. The amplitude of nonaxisymmetry and nonstationarity and its impact on the K_{z} studies is in general accepted to be small (Read 2014), but should be quantified precisely on an observational basis. Much larger samples with extremely accurate data on a much wider Galactic volume are expected from the Gaia mission, and will help in examining and solving all such questions.
At very low b latitudes, not considered here, a colour cut criterion has been applied (see Kordopatis et al. 2013).
Acknowledgments
Funding for RAVE has been provided by the Australian Astronomical Observatory; the LeibnizInstitut für Astrophysik Potsdam (AIP); the Australian National University; the Australian Research Council; the French National Research Agency; the German Research Foundation (SPP 1177 and SFB 881); the European Research Council (ERCStG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; Johns Hopkins University; the National Science Foundation of the USA (AST0908326); the W. M. Keck foundation; the Macquarie University; the Netherlands Research School for Astronomy; the Natural Sciences and Engineering Research Council of Canada; the Slovenian Research Agency; the Swiss National Science Foundation; the Science & Technology Facilities Council of the UK; Opticon; Strasbourg Observatory and the universities of Groningen, Heidelberg, and Sydney. The RAVE website is at http://www.ravesurvey.org. We thank Annie Robin and Michel Crézé for providing Fig. 11 and for comments.
References
 Abadi, M. G., Navarro, J. F., Fardal, M., Babul, A., & Steinmetz, M. 2010, MNRAS, 407, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Batsleer, P., & Dejonghe, H. 1994, A&A, 287, 43 [NASA ADS] [Google Scholar]
 Bienaymé, O. 1999, A&A, 341, 86 [NASA ADS] [Google Scholar]
 Bienaymé, O., Soubiran, C., Mishenina, T. V., Kovtyukh, V. V., & Siebert, A. 2006, A&A, 446, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bienaymé, O., Famaey, B., Wu, X., Zhao, H. S., & Aubert, D. 2009, A&A, 500, 801 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Binney, J. 2012, MNRAS, 426, 1324 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., Bissantz, N., & Gerhard, O. 2000, ApJ, 537, L99 [NASA ADS] [CrossRef] [Google Scholar]
 Bissantz, N., Englmaier, P., & Gerhard, O. 2003, MNRAS, 340, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Bobylev, V. V., & Bajkova, A. T. 2013, Astron. Lett., 39, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., & Rix, H.W. 2013, ApJ, 779, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Bovy, J., Allen de Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Brunthaler, A., Reid, M. J., Menten, K. M., et al. 2011, Astron. Nachr., 332, 461 [NASA ADS] [CrossRef] [Google Scholar]
 Büdenbender, A., van de Ven, G., & Watkins, L. L. 2014, MNRAS, submitted [arXiv:1407.4808] [Google Scholar]
 Burch, B., & Cowsik, R. 2013, ApJ, 779, 35 [NASA ADS] [CrossRef] [Google Scholar]
 CabreraLavers, A., Garzón, F., & Hammersley, P. L. 2005, A&A, 433, 173 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cannon, R. D., & Lloyd, C. 1969, MNRAS, 144, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Chang, C.K., Ko, C.M., & Peng, T.H. 2011, ApJ, 740, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Crézé, M., Chereul, E., Bienaymé, O., & Pichon, C. 1998, A&A, 329, 920 [NASA ADS] [Google Scholar]
 Debattista, V. P. 2014, MNRAS, 443, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Dejonghe, H., & de Zeeuw, T. 1988, ApJ, 333, 90 [NASA ADS] [CrossRef] [Google Scholar]
 de Zeeuw, T. 1985, MNRAS, 216, 273 [NASA ADS] [CrossRef] [Google Scholar]
 ESA 1997, The Hipparcos catalogue, ESA SP, 1200 [Google Scholar]
 Famaey, B., & Binney, J. 2005, MNRAS, 363, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & Dejonghe, H. 2003, MNRAS, 340, 752 [NASA ADS] [CrossRef] [Google Scholar]
 Famaey, B., & McGaugh, S. S. 2012, Liv. Rev. Rel., 15, 10 [Google Scholar]
 Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564 [NASA ADS] [CrossRef] [Google Scholar]
 Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [CrossRef] [Google Scholar]
 Garbari, S., Liu, C., Read, J. I., & Lake, G. 2012, MNRAS, 425, 1445 [NASA ADS] [CrossRef] [Google Scholar]
 Girardi, L. 1999, MNRAS, 308, 818 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Groenewegen, M. A. T. 2008, A&A, 488, 935 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hammer, F., Puech, M., Flores, H., et al. 2012, Eur. Phys. J. Web Conf., 19, 1004 [CrossRef] [EDP Sciences] [Google Scholar]
 Holmberg, J., & Flynn, C. 2004, MNRAS, 352, 440 [NASA ADS] [CrossRef] [Google Scholar]
 Holmberg, J., Flynn, C., & Portinari, L. 2006, MNRAS, 367, 449 [NASA ADS] [CrossRef] [Google Scholar]
 James, F. 2004, Reprinted from the Proceedings of the 1972 CERN Computing and Data Processing School [Google Scholar]
 Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Ibata, R., Lewis, G. F., Martin, N. F., Bellazzini, M., & Correnti, M. 2013, ApJ, 765, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Kapteyn, J. C. 1922, ApJ, 55, 302 [NASA ADS] [CrossRef] [Google Scholar]
 Korchagin, V. I., Girard, T. M., Borkova, T. V., Dinescu, D. I., & van Altena, W. F. 2003, AJ, 126, 2896 [NASA ADS] [CrossRef] [Google Scholar]
 Koposov, S. E., Rix, H.W., & Hogg, D. W. 2010, ApJ, 712, 260 [NASA ADS] [CrossRef] [Google Scholar]
 Kordopatis, G., Gilmore, G., Steinmetz, M., et al. 2013, AJ, 146, 134 [NASA ADS] [CrossRef] [Google Scholar]
 Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571 [NASA ADS] [CrossRef] [Google Scholar]
 Kuzmin, G. G., & Kutuzov, S. A. 1962, Bull. Abastumani Ap. Obs., 27, 82 [Google Scholar]
 Laney, C. D., Joner, M. D., & Pietrzyński, G. 2012, MNRAS, 419, 1637 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, J. R., & Freeman, K. C. 1989, AJ, 97, 139 [NASA ADS] [CrossRef] [Google Scholar]
 Ling, F.S. 2010, Phys. Rev. D, 82, 3534 [NASA ADS] [CrossRef] [Google Scholar]
 Macciò, A. V., Dutton, A. A., van den Bosch, F. C., et al. 2007, MNRAS, 378, 55 [NASA ADS] [CrossRef] [Google Scholar]
 Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940 [NASA ADS] [CrossRef] [Google Scholar]
 McMillan, P. J. 2011, MNRAS, 414, 2446 [NASA ADS] [CrossRef] [Google Scholar]
 Miglio, A., Chiappini, C., Morel, T., et al. 2013, MNRAS, 429, 423 [NASA ADS] [CrossRef] [Google Scholar]
 Mignard, F. 2000, A&A, 354, 522 [NASA ADS] [Google Scholar]
 Milgrom, M. 1983, ApJ, 270, 365 [NASA ADS] [CrossRef] [Google Scholar]
 Nataf, D. M., Gould, A., & Pinsonneault, M. H. 2012, Acta Astron., 62, 33 [NASA ADS] [Google Scholar]
 Navarro, J. F., & Steinmetz, M. 2000, ApJ, 528, 607 [NASA ADS] [CrossRef] [Google Scholar]
 Nesti, F., & Salucci, P. 2013, J. Cosmol. Astropart. Phys., 7, 16 [Google Scholar]
 Ojha, D. K. 2001, MNRAS, 322, 426 [NASA ADS] [CrossRef] [Google Scholar]
 Ojha, D. K., Bienaymé, O., Robin, A. C., Crézé, M., & Mohan, V. 1996, A&A, 311, 456 [NASA ADS] [Google Scholar]
 Oort, J. H. 1932, Bull. Astron. Inst. Netherlands, 6, 249 [Google Scholar]
 Olling, R. P., & Dehnen, W. 2003, ApJ, 599, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Olling, R. P., & Merrifield, M. R. 1998, MNRAS, 297, 943 [NASA ADS] [CrossRef] [Google Scholar]
 Pasetto, S., Grebel, E. K., Zwitter, T., et al. 2012a, A&A, 547, A71 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pasetto, S., Grebel, E. K., Zwitter, T., et al. 2012b, A&A, 547, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piffl, T., Scannapieco, C., Binney, J., et al. 2014a, A&A, 562, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Piffl, T., Binney, J., McMillan, P., et al. 2014b, MNRAS, submitted [arXiv:1406:4130v1] [Google Scholar]
 Pillepich, A., Kuhlen, M., Guedes, J., & Madau, P. 2014, ApJ, 784, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Polido, P., Jablonski, F., & Lépine, J. R. D. 2013, ApJ, 778, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Read, J. I. 2014, J. Phys. G Nucl. Phys., 41, 3101 [Google Scholar]
 Read, J. I., Mayer, L., Brooks, A. M., Governato, F., & Lake, G. 2009, MNRAS, 397, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130 [Google Scholar]
 Salucci, P., Nesti, F., Gentile, G., & Frigerio Martins, C. 2010, A&A, 523, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Siebert, A., Bienaymé, O., & Soubiran, C. 2003, A&A, 399, 531 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Siebert, A., Bienaymé, O., Binney, J., et al. 2008, MNRAS, 391, 793 [NASA ADS] [CrossRef] [Google Scholar]
 Siebert, A., Famaey, B., Minchev, I., et al. 2011a, MNRAS, 412, 2026 [NASA ADS] [CrossRef] [Google Scholar]
 Siebert, A., Williams, M. E. K., Siviero, A., et al. 2011b, AJ, 141, 187 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H. 1969, ApJ, 158, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, MNRAS, 379, 755 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, M. C., Whiteoak, S. H., & Evans, N. W. 2012, ApJ, 746, 181 [NASA ADS] [CrossRef] [Google Scholar]
 Soubiran, C., Bienaymé, O., & Siebert, A. 2003, A&A, 398, 141 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Soubiran, C., Bienaymé, O., Mishenina, T. V., & Kovtyukh, V. V. 2008, A&A, 480, 91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Statler, T. S. 1989, ApJ, 344, 217 [NASA ADS] [CrossRef] [Google Scholar]
 Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645 [NASA ADS] [CrossRef] [Google Scholar]
 Stello, D., Huber, D., Bedding, T. R., et al. 2013, ApJ, 765, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Valentini, M., & Munari, U. 2010, A&A, 522, A79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 van der Kruit, P. C., & Freeman, K. C. 2011, ARA&A, 49, 301 [NASA ADS] [CrossRef] [Google Scholar]
 van Leeuwen, F. 2007, A&A, 474, 653 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Veltz, L., Bienaymé, O., Freeman, K. C., et al. 2008, A&A, 480, 753 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, L., Rix, H.W., van de Ven, G., et al. 2013, ApJ, 772, 108 [NASA ADS] [CrossRef] [Google Scholar]
 Zhao, G., Qiu, H.M., & Zhang, H.W. 2000, Acta Astrophys. Sinica, 20, 389 [NASA ADS] [Google Scholar]
 Zwitter, T., Siebert, A., Munari, U., et al. 2008, AJ, 136, 421 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Separable potential and Jeans equation of the vertical motion
(Summary) Here we show that the Jeans equation of the vertical motion of stars and the Collisionless Boltzmann equation (CBE) for the vertical motions have the same general solutions. The CBE is obtained under the assumption of separability of the vertical and horizontal motions (i.e. σ_{r,z} = 0), while the Jeans equation is obtained under a less restrictive assumption.
Garbari et al. (2012) used a general formulation for the solutions of the Jeans equation of vertical motion. Then, they claim that their solution is more general than that obtained under the more restrictive hypothesis of separability.
Here, we show the opposite. Under these two different hypotheses (one more restrictive than the other), we obtain the same general solutions. Thus, the useful formulation used by Garbari et al. (2012) is not more general than that obtained under the simple hypothesis of separability.
(Justification) The usual approximation, at low z, to model the motion of stars perpendicular to the Galactic plane consists of assuming that the Galactic potential is separable in R and z coordinates, thus the correlation between vertical and horizontal velocities is zero, σ_{Rz} = 0 and the velocity ellipsoids remain parallel to the Galactic plane. Reciprocally, the potential is separable in the domains where σ_{Rz} = 0.
Under the assumption of separability, the stationary vertical distribution function of stars may be described using the 2D collisionless Boltzmann equation (CBE): $\mathit{w}\hspace{0.17em}\frac{\mathit{\partial f}\mathrm{\left(}\mathit{z,w}\mathrm{\right)}}{\mathit{\partial z}}\mathrm{}\frac{\mathit{\partial}\mathrm{\Phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{\partial z}}\hspace{0.17em}\frac{\mathit{\partial f}\mathrm{\left(}\mathit{z,w}\mathrm{\right)}}{\mathit{\partial w}}\mathrm{=}\mathrm{0.}$(A.1)Besides, the Jeans equation corresponding to stationary vertical motions is $\frac{\mathrm{1}}{\mathit{R}}\hspace{0.17em}\frac{\mathit{\partial}}{\mathit{\partial R}}\mathrm{(}\mathit{R\nu}{{\mathit{\sigma}}_{\mathit{Rz}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{(}\mathit{\nu}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathit{\nu}\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial z}}\mathrm{=}\mathrm{0}$(A.2)and can also be simplified by canceling the first l.h.s. term when the potential is separable since then σ_{Rz} = 0.
Recently, Garbari et al. (2012) used a general solution of this reduced Jeans equation: $\frac{\mathit{\partial}}{\mathit{\partial z}}\mathrm{(}\mathit{\nu}{{\mathit{\sigma}}_{\mathit{z}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathit{\nu}\frac{\mathit{\partial}\mathrm{\Phi}}{\mathit{\partial z}}\mathrm{=}\mathrm{0.}$(A.3)They notice that the cancellation of the first l.h.s. term of Eq. (A.2) covers more general and less restricting hypothesis than the assumption of separability (i.e. σ_{Rz} = 0) and claim that this “minimal assumption method ... breaks the assumption that the distribution is separable”. Without a clear justification, they assume that their general solution is more general than those usually used, for instance by Holmberg & Flynn (2004).
Here, we show the opposite and we establish that cancelling the l.h.s. term of the Jeans Eq. (A.1) gives the same solution that in the case of separability of the R and z motions.
Thus, let be f(z = 0,w) = g^{∗}(p = w^{2}/ 2), the velocity distribution of a stationary solution at z = 0. This function is odd and can be written as an integration over an infinite set of Gaussians: $\begin{array}{ccc}\mathit{f}\mathrm{(}\mathit{z}\mathrm{=}\mathrm{0}\mathit{,w}\mathrm{)}\mathrm{=}{\mathit{g}}^{\mathrm{\ast}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{a}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{}{\mathit{w}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}\mathrm{\right)}}\hspace{0.17em}\mathrm{d}\mathit{\sigma}\mathit{.}& & \end{array}$(A.4)This is equivalent to the Laplace transform, with β = 1 /σ^{2} and a(σ) dσ = g(β)dβ: ${\mathit{g}}^{\mathrm{\ast}}\mathrm{\left(}\mathit{p}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{g}\mathrm{\left(}\mathit{\beta}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{}\mathit{\beta p}}\hspace{0.17em}\mathrm{d}\mathit{\beta}\mathrm{=}\mathrm{\mathcal{L}}\hspace{0.17em}\mathrm{\left[}\mathit{g}\mathrm{\right(}\mathit{\beta}\mathrm{\left)}\mathrm{\right]}\mathit{.}$(A.5)If the integral exists, for instance when f(z = 0,w) is null for large values of w, then g(β) exists, it is unique and is given by the inverse Laplace transform $\mathit{g}\mathrm{\left(}\mathit{\beta}\mathrm{\right)}\mathrm{=}{\mathrm{\mathcal{L}}}^{1}\mathrm{\left[}{\mathit{g}}^{\mathrm{\ast}}\mathrm{\right(}\mathit{p}\mathrm{\left)}\mathrm{\right]}$(A.6)that gives us the unique decomposition in Gaussians (Eq. (A.4)).
For an isothermal component (i.e. Gaussian in velocities with a dispersion σ) the solution of the Jeans Eq. (A.3) is ν_{σ}(z) = e^{−βΦ(z)}, with Φ(0) = 0, and since, from Eq. (A.4), we have $\begin{array}{ccc}\mathit{\nu}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{a}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\sigma}& & \end{array}$then the general solution can be written: $\mathit{\nu}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{a}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{}\mathrm{\Phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{/}{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}\mathit{\sigma}\mathit{.}$(A.7)This general solution of Eq. (A.3) has a different form, but is identical to the general solution given by Eq. (8) of Garbari et al. (2012). A key point is that for any given odd function f(z = 0,w) representing the distribution function at z = 0, there is a unique function f′(z,w), solution of the Jeans Eq. (A.3), that satisfies f′(z = 0,w) = f(z = 0,w).
Now, on the other hand, the general solution of the CBE is h(E) where E is the energy. The distribution function at z = 0, is given by h(E(z = 0)) = h(w^{2}/ 2). As has been done previously, we can inverse the function h and rewrite the general solution. Thus, we obtain an equivalent form of the general solution: $\mathit{f}\mathrm{\left(}\mathit{z,w}\mathrm{\right)}\mathrm{=}{\mathrm{\int}}_{\mathrm{0}}^{\mathrm{\infty}}\mathit{a}\mathrm{\left(}\mathit{\sigma}\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{}\mathrm{\Phi}\mathrm{\left(}\mathit{z}\mathrm{\right)}\mathit{/}{\mathit{\sigma}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{1}}{\sqrt{\mathrm{2}\mathit{\pi}}\mathit{\sigma}}{\mathrm{e}}^{\mathrm{}{\mathit{w}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}\mathrm{\right)}}\mathrm{d}\mathit{\sigma ,}$(A.8)where a(σ) is related through an inverse Laplace transform to h(E = 0).
By integrating f(z,w) over the w velocities, we exactly recover the general solution (Eq. (A.7)) of the Jeans Eq. (A.3).
Thus, the general solution obtained for the Jeans Eq. (A.3) has the same form as the general solution of the CBE for a separable potential.
In conclusion, the hypothesis that the l.h.s. term in Eq. (A.2) is null, is indeed a more general hypothesis than the hypothesis of separability. However, contrary to the Garbari et al. (2012) claim, the general solution of Eq. (A.3) is neither more general nor different than the solution (Eq. (A.8)) obtained from the hypothesis of separability.
All Tables
All Figures
Fig. 1 Histograms of J − K colours for 12308 RAVE stars observed towards the South Galactic Pole (black continuous line). Red clump stars (red dashed), subgiants (blue dashdotted line), dwarfs (green dotted line) identified from their gravity. 

In the text 
Fig. 2 Histogram of gravities for 12308 RAVE stars observed towards the South Galactic Pole. 

In the text 
Fig. 3 2MASS star counts towards the South Galactic Pole versus the K_{c} magnitude (green circle) corrected for Galactic latitude (RAVE counts towards both poles are indicated by black crosses, and RAVE counts of red clump stars are denoted by red plus signs). 

In the text 
Fig. 4 Vertical number density distribution of red clump stars towards the South Galactic Pole (black symbols: 200 pc binning, red symbols 50 pc binning). 

In the text 
Fig. 5 Distribution of errors in the vertical Galactic velocity for stars with  z  < 2000 pc (continuous black line), and for stars with 1300 <  z  < 2000 pc (dotted red line). 

In the text 
Fig. 6 Vertical (black symbols) and radial (red symbols) velocity dispersions: σ_{Vz}, σ_{VR}. Mean vertical velocity V_{z} (black line). 

In the text 
Fig. 7 Vertical number density (top) and velocity dispersions (bottom), data (symbols with error bars) and model (lines) for the three RAVE metallicity samples and the NGPHipparcos sample: low metallicity [ M/H ] ≤ −0.35 (black), intermediate −0.35 < [ M/H ] ≤ −0.15 (red), high −0.15 < [ M/H ] (green), NGP (blue). Dotted lines are the extrapolated model not fitted to data. 

In the text 
Fig. 8 Vertical number density (top) and velocity dispersions (bottom) for the full RAVE RC sample (symbols with error bars) and model (continuous line). 

In the text 
Fig. 9 K_{z} force and 1σ error intervals. 

In the text 
Fig. 10 Vertical mass density of the stellar thin discs from the Besançon model (red line) fitted with an exponential of scale height h = 229 pc (black line). 

In the text 
Fig. 11 Total surface mass density (black continous line) $\mathrm{\Sigma}\mathrm{(}\mathit{<}\mathrm{}\mathit{z}\mathrm{\left}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int}}}_{\mathrm{}\mathit{z}}^{\mathit{z}}{\mathit{\rho}}_{\mathrm{tot}}\mathrm{d}\mathit{z}$ at R_{0} split in a DM spherical component (red dashed line) and in a baryonic double exponential disc (blue dashdotted line). Dotted lines are 1σ error intervals. 

In the text 
Fig. 12 Probability distribution function for the local dark matter density ρ_{DM}(z = 0). 

In the text 
Fig. 13 Probability distribution functions for the total local mass density, i.e. the Oort limit ρ_{tot}(z = 0) (continous line), and the baryonic mass density ρ_{bar}(z = 0) (dotted line). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.