# Weighing the local dark matter with RAVE red clump stars

###### Key Words.:

Galaxies: kinematics and dynamics## Abstract

Context: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 M pc, and (ii) at 2 kpc Mpc. Assuming the solar Galactic radius at kpc, we deduce the local dark matter density Mpc and the baryonic surface mass density Mpc. Our results are in agreement with previously published 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. 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 galactic-scale modifications of gravity.

Aims:

Methods:

Results:

Conclusions:

## 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 problem (where 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 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 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. 2011), 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. 2011; 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 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 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 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 problem and related questions by Read (2014).

This paper proceeds as follows. In Section 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 Section 4, we present the results. Finally, conclusions are given and discussed in Sections 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. Cross-identification 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 pc-radius 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 , and to minimize the Malmquist-type bias. To minimize the consequence of interstellar extinction, we also restrict the selection to Galactic directions . Thus, taking into consideration the extinction = 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 within the range from 0.5 to 0.8, based on 2MASS magnitudes (Figure 1).

Based on high-resolution 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 shows a well-defined peak, allowing us to distinguish RC stars from subgiants and dwarfs (Figure 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 figure 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 within 1.8 and 2.8, reduces the sample by 20%. This does not modify our final conclusions concerning the measurement of the force. It is a complementary indication that these stars, with RAVE are also RC stars.

Finally, we only select stars with a proper motion accuracy better than 4 mas y and better than 5 km s for the radial velocity, with S/N greater than 20 in the RAVE spectra, and . 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 with a dispersion of 0.22 magnitude, 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 . We note that the mean absolute -magnitude of RC stars does not depend on metallicity (Groenewegen 2008).

To define our main sample, we apply the supplementary selection criterion:

(1) |

This criterion includes RC stars located within the 500 pc-radius cylindric volume oriented towards the Galactic poles. The sample contains 9522 stars with 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 colour, apparent magnitudes (Eq. 1), and Galactic directions are applied.

RAVE is a magnitude-limited survey of stars randomly selected in the southern celestial hemisphere. The original design was to only observe stars in the interval , 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}

We determine this ratio using RAVE counts in both north and south Galactic directions in order to improve the statistics at low . The 2MASS counts are complete down to magnitude (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 , we do not determine the counts using the apparent magnitudes but using the corrected apparent magnitudes:

(2) |

which depends only on the absolute magnitude of the star and of its position.

Figure 3 plots the 2MASS counts as a function of magnitudes. Star counts from our main RAVE sample and from our RAVE sub-sample 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 , or equivalently as a function of the height (Figure 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 (), the accuracy of the various measured or used parameters remains high independent of the distance. For instance, the median accuracy in colours (within 0.5-0.8) is 0.03 from =6 to =10. Similarly, the mean S/N of the RAVE spectra used to determine the gravity remains high for RC stars at 2 kpc (10): the mean S/N is 51 (r.m.s 16). This implies that our selections and cuts remain homogeneous independent of the distance .

A second effect is the Malmquist bias: this depends on , 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, , with pc and =0.2, at z=1000 pc the bias on the estimated distances is using a cone for the counts and is 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 Sections 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, and , uncorrected for the solar motion, assuming =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% (Figure 5). The median error on the component is 2.4 km.s.

The mean vertical velocity is constant with (Figure 6). The velocity dispersions and are measured by applying a 3.5-sigma-clipping to the , Galactic velocity components. The uncertainties on the dispersions are . The vertical velocity dispersion rises up to 38 km s at 1 kpc and then remains nearly constant (Figure 6).

The velocity ellipsoid tilt is null at z=300 pc and reaches 81 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 =1kpc, it is about at 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 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 Section 3).

### 2.4 Metallicities

To improve the analysis of the vertical potential and the force, we split our sample according to the metallicity in three sub-samples 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 is mainly sensitive to the gradients of the density, and thus each sub-sample can probe more efficiently a different range of -distances, the lowest metallicity sample probing the potential for the largest distances.

On the other side, the vertical potential determination is poorly determined at low 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 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: = 2.08 for the metal-poor, low-mass end and reaches up to = 3 for the high-mass, metal-rich 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 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 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 century-old 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 mid-plane, the problem becomes 1D dynamical model. In this case, the vertical distribution of stellar positions and vertical velocities can be written as the sum of isothermal components depending only on the vertical energy:

with the vertical density of the tracer stars as:

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 , 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 Hamilton-Jacobi equation is fully separable, and a phase-space 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 and axis ratios of both components. These four free parameters allow us to constrain locally the scale-height 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 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 =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 of the foci along the vertical axis, defining a confocal ellipsoidal coordinate system. We set =2 kpc, in order to have a velocity ellipsoid at (, )=(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 Kuzmin-Kutuzov 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 () by:

(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 axis at .

A general Stäckel potential takes the form:

(4) |

where is an arbitrary function.

Here, we define the class of Kuzmin-Kutuzov potentials (Dejonghe & de Zeeuw 1988), writing (with the condition ):

(5) |

The corresponding isodensity surfaces from Poisson equation are flattened oblate spheroids. Increasing 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 remains the same.

Besides the energy and the angular momentum, a third independent isolating integral of the motion exists and can be written as:

(6) |

with

(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

(8) |

with the radius of the circular orbit that has the angular momentum , is the angular velocity, is the epicyclic frequency, and is the energy of a circular orbiting star at radius .

For sufficiently small velocity dispersions, the number density distribution, , is close to .

We set and the velocity dispersions are close to .

Local number density and dispersion, , , are constants, is (close to) the scale length of the number density distribution, and 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 . 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 =0, the corresponding density varies nearly exponentially in an extended domain of a few kpc around the Sun. We note that if 0, the density above the plane also depends on the vertical potential; thus the density may vary exponentially at any with a supplementary condition, as for instance (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: =2.5 kpc and =9 kpc. This is to be compared with 2.2 and 2.8 kpc for the thin and thick discs (Cabrera-Lavers et al. 2005; Jurić et al. 2008; Chang et al. 2011; Polido et al. 2013; Robin et al. 2014) and =4.4 or 5.6 kpc, (=8.8 or 11.1 kpc), respectively by Lewis & Freeman (1989) and Ojha et al. (1996). We also set =2 for the thin disc stars (km.s) 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 that increases the observed vertical velocity dispersion. The second effect is related to the vertical bending of the stellar orbits. For stars seen at =1 or 2 kpc, their mean Galactic radius , when they cross the Galactic plane, is larger than , a position where the stellar density is lower than at the radius . These two effects lead to an overestimate of in a 1D model. A third effect is also related to the bending of orbits and to the radial gradient of , lowering towards the pole. This effect leads to an underestimate of the 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 and vertical velocity dispersions 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 :

with

where the 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 expression is introduced to impose a potential with a nearly flat rotation curve. The factor does not need to be large and the contribution of the corresponding term to the 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 between =7 and 15 kpc.

We use two quasi-isothermal components (Eq. 8) to model the stellar sample distribution function, the density and dispersion for each metallicity sample.

We determine the 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 variance-covariance matrix. To compute the posterior probability distribution function (PDF), we consider the likelihood and apply a Markov chain Monte-Carlo using the Metropolis-Hastings Algorithm (Foreman-Mackey et al. 2013). From this, we determine the marginal PDFs.

Figure 7 shows the binned data and the best fit model to the density and dispersion versus for the three metallicity RAVE samples and for the Hipparcos-NGP sample. For information, we also show this best-fit model and the binned data for density and vertical velocity dispersions without metallicity splitting (Figure 8).

## 4 Results

### 4.1 The 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 century-old 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 -derivatives of this potential. It is known that the is an ill-conditioned 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 Kuzmin-Kutuzov 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 three-parameter modelling introduced by Kuijken & Gilmore (1989) to describe the total Galactic disc mass surface density. Our model goes beyond the plane-parallel approximation and allows us to describe the distribution function 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.

The resulting at from the fit of RAVE data (Figures 7) is shown in Figure 9. We give the results in terms of G for visibility and comparison with other studies, even though this should not be confused with the true surface density.

The obtained is less constrained at low where there are no RAVE data, but only the smaller sample of Hipparcos and Elodie stars. The force at 350 pc is found to be 44.2Mpc in agreement with the Korchagin et al. (2003) determination Mpc, based on a sample of 1500 K giants from Hipparcos observations, mainly first ascent giants rather than clump giants.

We also note that our determination near =0 allows us to deduce the Oort limit (z=0)=0.0911Mpc, intermediate value between the Crézé et al. (1998) (0.076Mpc) and Holmberg & Flynn (2004) (0.102Mpc) determinations. These two studies are based on the same Hipparcos data, below =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 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 quasi-totality of the ordinary matter resides below =1 kpc, the mass density beyond 1 kpc is dominated by the dark matter. This results that our measurement gives direct access to the dark matter density between 1 and 2 kpc above the Galactic plane.

The force at intermediate and high distances are:

(1kpc)/(2G) =68.51.0pc =185227 km s kpc, and

(2kpc)/(2G) =96.92.2pc =261959 km s kpc .

### 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 at =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, non-existent 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, , modifies by less than 1% the between = 0 to 2 kpc, when is varying from 1.8 kpc to 3.5 kpc. Increasing the kinematics scale length does not modify the determination. Only decreasing the scale length to 7 or 5 kpc (or to 3.5 or 2.5 kpc) increases the 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 determination is not completely sufficient, and we must also know the 3D shape of the baryonic and dark matter components. Here, the four-parameter Stäckel potential we fitted should be considered simply as a way to estimate the 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 , 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. 2014) 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 (Figure 10).

For a given mass surface density of the disc at , both parameters and can modify the vertical force. Decreasing the scale length increases the vertical force of the disc. In this case, to fit the observed , the surface mass density of the disc must be decreased (and to fit the at larger , the dark matter mass density must be increased).

Here, we will consider as reference values =2200 pc, =300 pc and =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 double-exponential 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 . The resulting total vertical mass density distribution is shown in Figure 11 with the DM halo and disc decomposition. The plotted surface densities are the integrated mass volume density between and . We obtain for the total surface density of the disc component at the solar position:

Mpc.

The probability distribution function of total local volume mass density of the dark matter component is plotted in Figure 12. This yields

=0)=0.01430.0011 Mpc = 0.540.004 .

In the case of the baryonic and the total (baryonic+DM) local volume densities, the Oort limit, we obtain

=Mpc,

=Mpc,

and Figure 13 plots their respective probability distribution function.

Table 1 shows the results for and in the case of some other disc parameter values (results are based on a minimization and are slightly different from these resulting from MCMCs). The change of the best-fit 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.

pc | pc | pc | km.s | ||
---|---|---|---|---|---|

8500 | 2200 | 300 | 45.6 | 0.0154 | 267 |

8500 | 2000 | 300 | 45.0 | 0.0157 | 263 |

8500 | 2400 | 300 | 46.5 | 0.0146 | 264 |

8500 | 2200 | 250 | 43.2 | 0.0158 | 268 |

7500 | 2200 | 300 | 66.4 | 0.0166 | 263 |

We can note that, remarkably, the dark matter density is nearly insensitive to the model parameters , , . 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 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 problem is, in principle, ill-conditioned. 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 three-parameter 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 determination, but allows us a more realistic dynamical self-consistent distribution function modelling the density and velocity distribution of tracer stars.

### 5.1 The force measurement

At =1 kpc, our measurement is similar and in agreement with the two last decades determinations of the force based on data from stars up to 1 kpc: kpc/2 G)68 M pc (for recent determinations, see Garbari et al. 2012; Bovy & Rix 2013; Zhang et al. 2013). Between =1 and 2 kpc, we find that the scale-height of the dark halo is significantly larger than the domain probed with our sample, and that the density is nearly constant below =2 kpc.

The local density of the halo is =0)= 0.01430.0011M pc=0.5420.042 Gev cm. We also determine the total surface density of the disc components =44.44.1 M pc 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 (z=0)=0.006-0.008 M pc (see the compilation in Table 4 by Read 2014) However, it must be noted that in these previous measurements, the 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 non-local Galactic characteristics as the scale length of the disc, the Galactic solar radius , 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 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 spectro-photometric distances (Bovy et al. 2012) that favour a flat rotation curve in an extended interval of radii around .

Now, considering such a flat rotation curve with 220 km s, 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 0.006-0.01 M pc 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 240 km s(McMillan 2011; Bobylev & Bajkova 2013; Reid et al. 2014), and also assuming a spherical dark matter halo, McMillan (2011) obtained =0.400.04 Gev cm(=0.010 M/pc) 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 = km s kpc (and =+0.02). Their value of the slope of velocity curve introduces a small additional non-zero contribution to the local mass density from the Poisson equation:

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 at the Galactic centre and at from published observations of the circular velocity curve. They plot the force that can be directly compared with our measurement. One of their models (figure 6b) is in agreement with our determination within =0 to 2 kpc and they obtained Mpc. Unfortunately, in their model they adopt a rapidly rising rotation curve at that we judge very unlikely. This rising rotation curve, with 0.16 explains their high value obtained for . In their model, if the rotation curve was flat, the local density would probably be about 0.006 Mpc.

In a recent study, using tracer stars with between 1 to 2 kpc, Smith et al. (2012) estimated an order of magnitude for M pc, 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 .

### 5.3 Global Galactic properties

When we consider our model with =0.014 M pc, a flat rotation curve and a spherical halo, this implies =267 km s, 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 =220 km s 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 =0.99, instead of =0.53 in the case of DM-only numerical simulations.

In fact, combining the circular velocity from the BeSSeL project 240 km s(McMillan 2011; Bobylev & Bajkova 2013; Reid et al. 2014), with a nearly flat rotation curve at the Sun position, and a flattening of , 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 Tully-Fisher 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, figure 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 non-rotating 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 force probably does not include a sufficient number of free parameters to accurately model the shape of the force between 1 and 2 kpc. We might suspect that our modelling forces a nearly linear rise of the between 1 and 2 kpc. A supplementary (thick disc) Kuzmin-Kutuzov 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 /(2G) 90 M pc is predicted at =2 kpc. Nevertheless, a somewhat high value 75 M pc is rather predicted at kpc. For the same reason mentioned previously, our two Kuzmin-Kotozov components modelling of the vertical force is not adequate to correctly reproduce the force in the case of Mondian model that shows a significant bending of the vertical force law between =1 and 2 kpc.

### 5.4 A Galactic DM halo with core?

Now, if we consider the consistency of our local 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 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 (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 . This has implications for the mass concentration and can be tested relative to the mass-concentration 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 (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 =0)= 0.01430.0011M pc. We have independently determined the local DM density and the baryonic disc surface density at the solar Galactic radius . 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 beyond 1 kpc and are not modelled by Stäckel potentials, or with a simple torus fitting. Also, the non-axisymmetric effects due to spiral arms (Faure et al. 2014; Debattista 2014) or non-equilibrium 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 non-axisymmetry and non-stationarity and its impact on the 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.

###### Acknowledgements.

Funding for RAVE has been provided by the Australian Astronomical Observatory; the Leibniz-Institut 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 (ERC-StG 240271 Galactica); the Istituto Nazionale di Astrofisica at Padova; Johns Hopkins University; the National Science Foundation of the USA (AST-0908326); 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.rave-survey.org. We thank Annie Robin and Michel Crézé for providing figure 11 and for comments.## 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. ), 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 , to model the motion of stars perpendicular to the Galactic plane consists of assuming that the Galactic potential is separable in and coordinates, thus the correlation between vertical and horizontal velocities is zero, and the velocity ellipsoids remain parallel to the Galactic plane. Reciprocally, the potential is separable in the domains where .

Under the assumption of separability, the stationary vertical distribution function of stars may be described using the 2D collisionless Boltzmann equation (CBE):

(9) |

Besides, the Jeans equation corresponding to stationary vertical motions is

(10) |

and can also be simplified by canceling the first l.h.s. term when the potential is separable since then .

Recently, Garbari et al. (2012) used a general solution of this reduced Jeans equation:

(11) |

They notice that the cancellation of the first l.h.s. term of Equation 10 covers more general and less restricting hypothesis than the assumption of separability (i.e. ) 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 Equation 9 gives the same solution that
in the case of separability of the and motions.

Thus, let be =0,) = =, the velocity distribution of a stationary solution at =0. This function is odd and can be written as an integration over an infinite set of Gaussians:

(12) |

This is equivalent to the Laplace transform, with and :

(13) |

If the integral exists, for instance when is null for large values of , then exists, it is unique and is given by the inverse Laplace transform

(14) |

that gives us the unique decomposition in Gaussians (Equation 12).

For an isothermal component (i.e. Gaussian in velocities with a dispersion ) the solution of the Jeans Equation 11 is , with (0)=0, and since, from Equation 12, we have

then the general solution can be written:

(15) |

This general solution of Equation 11 has a different form, but is identical to the general solution given by the equation 8 of Garbari et al. (2012). A key point is that for any given odd function =0, representing the distribution function at =0,
there is a unique function , solution of the Jeans Equation 11, that satisfies =0,) = =.

Now, on the other hand, the general solution of the CBE is where is the energy. The distribution function at =0, is given by =0)) = . As has been done previously, we can inverse the function and rewrite the general solution. Thus, we obtain an equivalent form of the general solution:

(16) |

where is related through an inverse Laplace transform to =0).

By integrating over the velocities, we exactly recover the general solution (Equation 15) of the Jeans Equation 11.

Thus, the general solution obtained for the Jeans Equation 11 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 Equation 10 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 Equation 11 is neither more general nor different than the solution (Equation 16) obtained from the hypothesis of separability.

### Footnotes

- At very low latitudes, not considered here, a colour cut criterion has been applied (see Kordopatis et al. 2013)

### References

- Abadi, M. G., Navarro, J. F., Fardal, M., Babul, A., & Steinmetz, M. 2010, MNRAS, 407, 435
- Batsleer, P., & Dejonghe, H. 1994, A&A, 287, 43
- Bienaymé, O. 1999, A&A, 341, 86
- Bienaymé, O., Famaey, B., Wu, X., Zhao, H. S., & Aubert, D. 2009, A&A, 500, 801
- Bienaymé, O., Soubiran, C., Mishenina, T. V., Kovtyukh, V. V., & Siebert, A. 2006, A&A, 446, 933
- Binney, J., Bissantz, N., & Gerhard, O. 2000, ApJ, 537, L99
- Binney, J. 2012, MNRAS, 426, 1324
- Bissantz, N., Englmaier, P., & Gerhard, O. 2003, MNRAS, 340, 949
- Bobylev, V. V., & Bajkova, A. T. 2013, Astronomy Letters, 39, 809
- Bovy, J., & Rix, H.-W. 2013, ApJ, 779, 115
- Bovy, J., Allende Prieto, C., Beers, T. C., et al. 2012, ApJ, 759, 131
- Brunthaler, A., Reid, M. J., Menten, K. M., et al. 2011, Astronomische Nachrichten, 332, 461
- Büdenbender, A., van de Ven, G., & Watkins, L. L. 2014, arXiv:1407.4808
- Burch, B., & Cowsik, R. 2013, ApJ, 779, 35
- Cabrera-Lavers, A., Garzón, F., & Hammersley, P. L. 2005, A&A, 433, 173
- Cannon, R. D., & Lloyd, C. 1969, MNRAS, 144, 449
- Chang, C.-K., Ko, C.-M., & Peng, T.-H. 2011, ApJ, 740, 34
- Crézé, M., Chereul, E., Bienaymé, O., & Pichon, C. 1998, A&A, 329, 920
- Debattista, V. P. 2014, arXiv:1405.6345
- Dejonghe, H., & de Zeeuw, T. 1988, ApJ, 333, 90
- de Zeeuw, T. 1985, MNRAS, 216, 273
- ESA 1997, The Hipparcos catalogue, ESA SP-1200
- Famaey, B., & Dejonghe, H. 2003, MNRAS, 340, 752
- Famaey, B., & Binney, J. 2005, MNRAS, 363, 603
- Famaey, B., & McGaugh, S. S. 2012, Living Reviews in Relativity, 15, 10
- Faure, C., Siebert, A., & Famaey, B. 2014, MNRAS, 440, 2564
- Flynn, C., Holmberg, J., Portinari, L., Fuchs, B., & Jahreiß, H. 2006, MNRAS, 372, 1149
- Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
- Garbari, S., Liu, C., Read, J. I., & Lake, G. 2012, MNRAS, 425, 1445
- Girardi, L. 1999, MNRAS, 308, 818
- Girardi, L., Bressan, A., Bertelli, G., & Chiosi, C. 2000, A&AS, 141, 371
- Gómez, F. A., Minchev, I., O’Shea, B. W., et al. 2013, MNRAS, 429, 159
- Groenewegen, M. A. T. 2008, A&A, 488, 935
- Hammer, F., Puech, M., Flores, H., et al. 2012, European Physical Journal Web of Conferences, 19, 1004
- Holmberg, J., & Flynn, C. 2004, MNRAS, 352, 440
- Holmberg, J., Flynn, C., & Portinari, L. 2006, MNRAS, 367, 449
- James, F. 2004, Reprinted from the Proceedings of the 1972 CERN Computing and Data Processing School
- Jurić, M., Ivezić, Ž., Brooks, A., et al. 2008, ApJ, 673, 864
- Ibata, R., Lewis, G. F., Martin, N. F., Bellazzini, M., & Correnti, M. 2013, ApJ, 765, L15
- Kapteyn, J. C. 1922, ApJ, 55, 302
- Korchagin, V. I., Girard, T. M., Borkova, T. V., Dinescu, D. I., & van Altena, W. F. 2003, AJ, 126, 2896
- Koposov, S. E., Rix, H.-W., & Hogg, D. W. 2010, ApJ, 712, 260
- Kuijken, K., & Gilmore, G. 1989, MNRAS, 239, 571
- Kuzmin, G. G., & Kutuzov, S. A. 1962, Bull Abastumani Ap. Obs., 27, 82
- Kordopatis, G., Gilmore, G., Steinmetz, M., et al. 2013, AJ, 146, 134
- Laney, C. D., Joner, M. D., & Pietrzyński, G. 2012, MNRAS, 419, 1637
- Lewis, J. R., & Freeman, K. C. 1989, AJ, 97, 139
- Ling, F.-S. 2010, Phys. Rev. D, 82, 023534
- Macciò, A. V., Dutton, A. A., van den Bosch, F. C., et al. 2007, MNRAS, 378, 55
- Macciò, A. V., Dutton, A. A., & van den Bosch, F. C. 2008, MNRAS, 391, 1940
- McMillan, P. J. 2011, MNRAS, 414, 2446
- Miglio, A., Chiappini, C., Morel, T., et al. 2013, MNRAS, 429, 423
- Mignard, F. 2000, A&A, 354, 522
- Milgrom, M. 1983, ApJ, 270, 365
- Nataf, D. M., Gould, A., & Pinsonneault, M. H. 2012, Acta Astron., 62, 33
- Navarro, J. F., & Steinmetz, M. 2000, ApJ, 528, 607
- Nesti, F., & Salucci, P. 2013, J. Cosmology Astropart. Phys., 7, 16
- Ojha, D. K. 2001, MNRAS, 322, 426
- Ojha, D. K., Bienaymé, O., Robin, A. C., Crézé, M., & Mohan, V. 1996, A&A, 311, 456
- Oort, J. H. 1932, Bull. Astron. Inst. Netherlands, 6, 249
- Olling, R. P., & Dehnen, W. 2003, ApJ, 599, 275
- Olling, R. P., & Merrifield, M. R. 1998, MNRAS, 297, 943
- Pasetto, S., Grebel, E. K., Zwitter, T., et al. 2012, A&A, 547, A71
- Pasetto, S., Grebel, E. K., Zwitter, T., et al. 2012, A&A, 547, A70
- Perryman, M. A. C., de Boer, K. S., Gilmore, G., et al. 2001, A&A, 369, 339
- Piffl, T., Scannapieco, C., Binney, J., et al. 2014a, A&A, 562, A91
- Piffl, T., Binney, J., McMillan, P., et al. 2014b, arXiv 1404:4130v1
- Pillepich, A., Kuhlen, M., Guedes, J., & Madau, P. 2014, ApJ, 784, 161
- Polido, P., Jablonski, F., & Lépine, J. R. D. 2013, ApJ, 778, 32
- Read, J. I. 2014, Journal of Physics G Nuclear Physics, 41, 063101
- Read, J. I., Mayer, L., Brooks, A. M., Governato, F., & Lake, G. 2009, MNRAS, 397, 44
- Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130
- Robin, A.C., 2014, in prep
- Salucci, P., Nesti, F., Gentile, G., & Frigerio Martins, C. 2010, A&A, 523, A83
- Siebert, A., Bienaymé, O., & Soubiran, C. 2003, A&A, 399, 531
- Siebert, A., Bienaymé, O., Binney, J., et al. 2008, MNRAS, 391, 793
- Siebert, A., Famaey, B., Minchev, I., et al. 2011, MNRAS, 412, 2026
- Siebert, A., Williams, M. E. K., Siviero, A., et al. 2011, AJ, 141, 187
- Shu, F. H. 1969, ApJ, 158, 505
- Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
- Smith, M. C., Ruchti, G. R., Helmi, A., et al. 2007, MNRAS, 379, 755
- Smith, M. C., Whiteoak, S. H., & Evans, N. W. 2012, ApJ, 746, 181
- Soubiran, C., Bienaymé, O., & Siebert, A. 2003, A&A, 398, 141
- Soubiran, C., Bienaymé, O., Mishenina, T. V., & Kovtyukh, V. V. 2008, A&A, 480, 91
- Statler, T. S. 1989, ApJ, 344, 217
- Steinmetz, M., Zwitter, T., Siebert, A., et al. 2006, AJ, 132, 1645
- Stello, D., Huber, D., Bedding, T. R., et al. 2013, ApJ, 765, L41
- Valentini, M., & Munari, U. 2010, A&A, 522, A79
- van der Kruit, P. C., & Freeman, K. C. 2011, ARA&A, 49, 301
- van Leeuwen, F. 2007, A&A, 474, 653
- Veltz, L., Bienaymé, O., Freeman, K. C., et al. 2008, A&A, 480, 753
- Williams, M. E. K., Steinmetz, M., Binney, J., et al. 2013, MNRAS, 436, 101
- Zhang, L., Rix, H.-W., van de Ven, G., et al. 2013, ApJ, 772, 108
- Zhao, G., Qiu, H.-M., & Zhang, H.-W. 2000, Acta Astrophysica Sinica, 20, 389
- Zwitter, T., Siebert, A., Munari, U., et al. 2008, AJ, 136, 421