Free Access
Issue
A&A
Volume 579, July 2015
Article Number A123
Number of page(s) 11
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/201526022
Published online 15 July 2015

© ESO, 2015

1. Introduction

Although there are many indirect lines of evidence for the presence of non-baryonic matter in the Universe, nowhere is a concrete estimate for the density of such dark matter (DM) more important than in the Solar neighbourhood: the local density determines how easily such matter might be detected in laboratory experiments. Starting with the classic analyses of Kapteyn (1922) and Oort (1932, 1960), there have been many attempts to use the kinematics of stars in the local Milky Way (MW) to estimate the local density and distribution of mass. Recently, the data required for such an analysis have become much better thanks to massive kinematic surveys like SEGUE (Yanny et al. 2009), RAVE (Steinmetz et al. 2006), and APOGEE (Allende Prieto et al. 2008). Read (2014) has recently given a detailed review of the theoretical expectations from a cosmological/extragalactic perspective and the history of the astronomical attempts to measure the local DM mass density; he shows that recent estimates have yielded values in the range ρDM = 5−15 mM pc-3, but that systematic effects are still the dominant source of error1. The measurement can be made using one of two different basic approaches, each with its own strengths and weaknesses.

One can attempt to fit the global kinematics of the stars and gas in the MW in detail and, by generating a global mass model, also estimate the local DM density (e.g. Bienaymé et al. 1987; McMillan 2011). The problem with this approach is that basic properties like the distances to kinematic structures (e.g. Binney & Merrifield 1998, Sect. 9.1.1; Reid et al. 2014) or even the distance to the center of the MW (e.g. Gillessen et al. 2009; Sofue et al. 2011; Branham 2014) are difficult to measure; the area within which a reasonably robust analysis can be made is only between the end of the central bar and the Sun’s Galactocentric radius, R0 (Dehnen & Binney 1998); and the parameters of arbitrary models for the DM geometry must be fit, e.g. using a Navarrow-Frenck-White (1997) or Einasto (1965) profile plus assumptions about the halo’s oblateness. A recent example of this approach is that of Piffl et al. (2014): using ~200 000 stars from the RAVE survey with Galactic heights | z | < 1.5 kpc and a 33-parameter model for the global mass-distribution and its effects on the stellar kinematics via assumed stellar distribution functions (DF), they derived a local DM density of 15 ± 2 mM pc-3. When kinematic data from the gaia satellite become available, this method will be the only way of dealing with the complexity of a global data sample from a realistic MW (Binney 1998) and the sheer enormity of the dataset will constrain the otherwise arbitrary model parameters to conform to reality.

An alternate approach better suited to the limits of current data and our knowledge of the structure of the MW (Carraro 2015) is to try to measure the DM density locally using local stellar kinematics (e.g. Bahcall 1984; Kuijken & Gilmore 1989c; Holmberg & Flynn 2004; Bovy & Rix 2013). The advantage of this approach is that the proper motions and distances to all stars and hence the full three-dimensional densities, velocities, and dispersions can principally be measured, at least at large distances from the plane of the MW where the effects of DM should show up most directly, and the large-scale distribution of DM presumedly plays a minor role (e.g. Siebert et al. 2008). Here too the complex potential and DF models can be used, but they may not be able to ensure that the result is accurate (versus merely precise). For example, the simple two-component Stäckel potentials used by Bienyamé et al. (2014) to model the disc and DM halo give the local disc a power-law rather than an exponential radial profile; any inaccuracies in the resulting model dispersion will be compensated by modifying the ill-constrained properties of the DM halo.

It is important to try to estimate the gravitational effects attributed to ρDM locally on length-scales smaller than R0, since there are competing theories like MOND (e.g. Milgrom 2014) that invoke changes in the behaviour of gravity on scales of whole galaxies – the same scales relevant for studies of the MW’s rotation curve. If the local effects of DM in the MW can only be determined using the MW’s rotation curve, we will have lost the information at one scale-length and have to deal with much more complexity and intrinsic uncertainties. The purpose of this paper is then to probe how robust truly local estimates for ρDM can be. I first show that the extent to which the derived values of ρDM are compromised by the many uncertainties and pitfalls in the methods and data used is generally not sufficently recognized. Then, I describe a relatively simple method with which one can principally make more robust estimates for the local DM density using the full vertical Jeans equation; as few assumptions as possible are made about the large-scale stellar kinematics, MW properties, and unknown DM model parameters, and the effects of radial coupling and a realistic rather than an idealized toy model is used for the local gravity. Finally, I apply this method to a particular dataset in order to see whether the challenges fundamentally lie in our ability to model the physics or to have adequate data. As we shall see, the problems with many previous ρDM estimates are serious and one must more carefully consider the inadequacies of models and the difficulty of defining and handling appropriate tracer datasets.

2. Near-midplane estimates for the DM density

If the stellar motions and the local potential Φ(z) near the MW mid-plane are simple enough, the Jeans theorem says that is an approximate isolating integral and the dynamic local density ρdyn(0) can be calculated from the distribution of mid-plane vertical velocities vz(0) and the observed density distribution ν(z) of a tracer population (Kuijken & Gilmore 1989c; see Sect. 4.9.3 in Binney & Tremaine 2008, hereafter BT08). Given enough good data, this method should yield a very robust estimate for ρdyn(0) and – given a sufficiently accurate model for the baryonic contributions – a similarly robust estimate for the remaining ρDM. This analysis has been made using many datasets and models (e.g. Kuijken & Gilmore 1989b; Holmberg & Flynn 2000; Creze et al. 1998; Garbari et al. 2012).

The problem with the conversion of the dynamic density, once measured, to the local DM density is the difficulty of estimating the baryonic contribution ρdyn,b(0). The observed mid-plane density of stars is fairly simple to estimate from local surveys: Read’s (2014) re-tabulation of the Flynn et al. (2006) model has 37.9 and 3.5 mM pc-3 from the thin and thick discs, respectively (if one assigns the local giants, white dwarfs, and brown dwarfs to the thin disc for this purpose). If one assumes that the density-ratio of the thick disc to thin disc is the value measured photometrically by Jurić et al. (2008), then the contributions have to be redistributed as 37.0 and 4.4 mM pc-3. The local surface density of molecular, neutral, and warm gas has been classically estimated as 13 M pc-2 by Holmberg & Flynn (2000), from which one can estimate the local densities using typical observed vertical scale-heights and/or linewidths. The stellar densities are generally thought to be good to about 10%, whereas the gas densities are difficult to estimate; the latter are given errors of ~50% by Homberg & Flynn (2000), suggesting a total error of about 25 mM pc-3. It should be noted that the local baryonic mass-density within this traditional model is dominated by the gas – i.e. by the least well-determined component – and the uncertainty is already much larger than the remaining 5−15 mM pc-3 thought to be due to DM.

The total estimated baryonic density ρb(0) ≈ 91 mM pc-3 can be converted to total surface density by adding up the contributions of individual kinematically isothermal components and hence depends upon the assumed or fitted potential. Using the potential from Holmberg & Flynn (2004), one obtains the values shown in Table 1. Like Read, I have used updated HI dispersions (Kalberla & Dedes 2008). These surface density estimates are not internally-consistent – one would have to repeat the fit of the K giant stars performed by Holmberg & Flynn using the modified densities and dispersions – but they are similar to those tabulated by Flynn et al. or Read and are good enough for our present purposes.

Table 1

A local baryonic mass model using classic (updated) ISM density estimates.

The traditional estimates of the neutral and molecular gaseous densities from the 1980 s and 1990 s are off for a variety of reasons; we now know that the density and ionization structure of the real ISM is much more complex. For example, a significant fraction of the HI cold neutral medium (CNM) is not optically thin, increasing the true densities by at least ~30% (Grenier et al. 2005) and probably by as much as factors of 23 (Braun 2012; Fukui et al. 2014). There is “dark” molecular gas not easily seen in CO but visible in gamma rays (e.g. Ackermann et al. 2012) and [CII] (e.g. Pineda et al. 2013) representing as much as 40% of the total, i.e. an increase from the traditional value of 67%. These effects must raise the estimates for the local ISM densities considerably: assuming a factor of 2 correction for the CNM and a factor of 1.67 for the molecular gas, ρgas,local(0) increases from 50 to 80 mM pc-3. The Sun also sits near the middle of the Local Bubble (e.g. Lallemont et al. 2003), a significant hole in the local Galactic interstellar medium (ISM) with a size comparable to the local thin-disc scale-height; the local gas density on somewhat larger scales is much higher by definition.

Unfortunately, the derived local value of ρb(0) is not even what we really need: ρb,dyn(0) and ρb,local(0) are not the same, given that the potential of the MW disc experienced by the stars is highly non-uniform and that the vertical epicylic frequencies are ~23 times larger than the orbital frequency (e.g. BT08, p. 167), effectively averaging out the local potential differences. The Sun sits in a trough between two major spiral arms (Perseus and Scutum-Centaurus) and at the edge of the Orion Spur (Xu et al. 2013): the local volume- and surface-densities must then be uncharacteristically lower than the mean kinematical values – again by definition – unless the effects of the Orion Spur are uncharacteristically large. It is extremely difficult to measure empirically or estimate theoretically how big this effect is, so it is generally not taken into account in the literature, not even when fitting local or global DM models. Drimmel & Spergel (2001) tried to measure the spiral amplitude of the MW by modelling projected COBE NIR data: they placed the Sun in an inter-arm region with a surface density contrast Σinterarm/ Σarm ≈ 0.85 in the K-band, corresponding to a spiral amplitude A = 0.14. External galaxies may give us a more reliable estimate for the typical arm-interarm contrast of large spiral galaxies like the MW (Rix & Zaritsky 1994). Grosbøl et al. (2004) analysed ground-based K-band images of 54 normal spiral galaxies and found amplitudes between 0.12 and 0.23. Elmegreen et al. (2011) and Kendall et al. (2015) studied 62 MW-like galaxies at 3.6 μm using the Spitzer satellite (Sheth et al. 2010; Kennicutt et al. 2003) and found typically larger ratios of A = 0.1−0.5. Thus, a galaxy like the MW is likely to have A ≈ 0.2, representing an interarm-to-mean surface density correction of 1/(1−A) ≈ 1.25. Vertical compression make the effects in volume density even larger than the ones in surface density: Dremmel & Spergel (2001) invoke a mid-plane arm-interarm flux density contrast that is 1.3 times larger than the projected value, yielding an interarm-to-mean volume density correction of 1/(1−1.3A) ≈ 1.35.

The effects of spiral arms are even stronger in the ISM: while Holwerda et al. (2005) find that there is only a difference of 2 in the opacities of arm and inter-arm regions, corresponding to the A ≈ 0.2 seen in external stellar discs, they attribute this effect to cloud clumpiness and not to total surface density. Indeed, the HI spiral amplitudes seen in the THINGS datasets (Walter et al. 2008) are much larger than those seen for stars, reaching up to nearly 100% amplitudes. Thus, it is not surprising that the mean HI surface densities at the Galactic radius of the Sun estimated by Kalberla & Dedes (2008) are larger than those of the local estimates2: (uncorrected for optical depth effects), consistent with a spiral amplitude A ≈ (15−11)/15 = 0.3.

Assuming that , a spiral amplitude of 20% alone increases the effective value of ρb,dyn(0) from 91 to 123 mM pc-3 and with more modern estimates of the densities to 163 mM pc-3. In order to compare these values of ρdyn(0) with those derived from recent kinematical studies (e.g. the exhaustive compilation in Read 2014), one can calculate the effective value of ρDM,effρdyn(0)−ρb,min(0) assuming Read’s baryonic minimum estimate of 91 mM pc-3: either the ISM and spiral amplitude corrections applied alone correspond to an apparent DM density of ~30 mM pc-3, and both corrections applied simultaneously result in ~72 mM pc-3. Assuming that the estimates for ρdyn(0) are accurate, it would appear that any reasonable correction to ρb(0) for known effects is much more than is needed and corresponds to apparent DM densities that are much larger than the range 5−15 mM pc-3 usually derived and thought to be consistent with cosmological expectations and galaxy-evolution scenarios. This is already true for the traditional estimate of the local baryonic density, given the uncertainties in the traditional estimate of the gas contribution, and the realistic corrections to the total ISM densities and/or corrections for spiral structure only make the situation much worse. Thus, practically any empirical value of ρdyn(0) in the range found by current kinematic studies is incapable of constraining the local DM density without additional (presently unavailable) information, assumptions, and/or additional external constraints.

3. Inhomogenous data and the effects of model assumptions

Since local estimates of ρDM using current data are unable to yield robust values by themselves, we must turn to kinematical effects farther away from the plane of the MW (e.g. Bahcall 1984; Garbari et al. 2011). By far the simplest method is to use the Jeans equations (moments of the Boltzmann equation; e.g. Sect. 4.8 in BT08), which connect the observed dispersions – particularly the zz component σzz(z) ≡ σz(z)2 – with the local gravity above and below the plane of the Galaxy. The latter is dominated by the relatively well-constrained stellar populations, small amounts of gas, and (presumedly) a roughly constant DM density. Here, there are also various approaches using different assumptions and amounts of data (e.g. Kuijken & Gilmore 1989a; Bovy & Tremaine 2012, hereafter BT12).

Moni Bidin et al. (2012b) analysed the kinematic data for a sample of colour-selected giants taken from Moni Bidin et al. (2012a; hereafter MBCM), for which the complete dispersion tensor could be estimated, using the Jeans and Poisson equations to derive the surface density with height, Σ(z). Assuming a double exponential tracer density distribution and a mean azimuthal velocity that is not a function of Galactocentric distance at any height from the plane, they derived a local DM density from the asymptotic behaviour of Σ(z) of 0 ± 1 mM pc-3. Bovy & Tremaine (2012) argued that one should assume a particular form for the vertical variation in the radial force field and a different radial scale-length, and so derived a value of 8 ± 3 mM pc-3 from a simple vertical Jeans equation (but see Moni Bidin et al. 2015). Beyond difficult technical questions about the measurability of dynamical quantities and their gradients and – in the case of the full Jeans equation – the values of various global MW parameters, all of these analyses suffer from a major defect: they incorrectly assume that the MBCM data represent the kinematics and density of a uniform tracer population with a well-defined uniform scale-length and scale-height. In fact, the giants are a vertically differentiated mixture of populations, each with its own very different abundance, scale-length, scale-height, and dispersion, and selection by colour alone is not enough to insure homogeneity of the spatially resolved kinematics. Since the vertical dispersion is a strong function of scale-height, the varying mixture of tracer populations automatically creates a large positive gradient in the mean σz as well as a decrease in the effective scale-height over the thick-disc value. Because the Jeans equations are linear, spatial and kinematic inhomogeneity is not per se a problem as long as the additional gradient in kinematics created by the mixture of tracer populations is accompanied by an equivalent change in the combined density. Thus, the impact of kinematic inhomogeneity depends upon the model used.

Equation (25) in BT12 – derived assuming strictly exponential density distributions with scale-heights h – can be used to understand the effects in the MBCM data at heights | z | ≫ h: assuming the change in the derived value of Σ(z) is mainly due to DM, that the effects of cross-dispersion terms are small, and that consistent with the current quality of the kinematic data – any gradient in is roughly constant, one can see that the DM reveals itself mostly via an induced vertical gradient in at large heights: (1)However, if this equation is applied to an inhomogeneous mixture of tracer populations rather than a homogeneous thick-disc tracer, an increase in the gradient of σz(z) appears due to the addition of tracers with lower σz and smaller h. If the density of this mixture is fit with a single exponential vertical density profile, the result will be a shorter apparent scale-height and the fitted DM density will be higher because of the increase in the apparent dispersion gradient and because of the decrease in apparent scale-height. One can easily estimate the magnitude of this effect using the simple thin+thick disc photometric model of Jurić et al. (2008): for two exponential components with scale-heights of 300 and 900 pc, a local thick-to-thin density ratio of 0.12, reasonable isothermal dispersions of 15 and 45 km s-1 (Kordopatis et al. 2013), and an observed value simply weighted by the relative densities, the tracer mixture effect alone is expected to produce a 20% increase in σz over the range of 1.4 to 4.5 kpc covered by the MBCM data, exactly as is observed. The corresponding change in effective scale-height depends upon the details of the sampling and fitting, but simple mock datasets created using the above model show a significant decrease in the effective vertical scale-height from the thick-disc value assumed by MB12. Bovy & Tremaine (2012) attempted to correct for the uncertainty in the mean scale-height, but if an inhomogeneous distribution is non-exponential, the fundamental equation and not the parameter is the problem. Sanders (2012) has performed a much more detailed simulation of the MBCM dataset and analysis using DF, a simple thin and thick disc, an assumed evolutionary history, a larger (20%) local thick disc fraction, and a standard DM halo. He – of course – found the same increase in the gradient in , but the complexity of his model has led to the impression that the increased gradient implies more rather than less DM (Garbari et al. 2012; Piffl et al. 2014).

A similar effect plagues the analysis of Bienaymé et al. (2014), who used RAVE observations of red clump giants, a simple potential model, and an assumed DF to derive the vertical gravity gz(z). Fitting gz(z) a posterori with a single doubly exponential disc, an intermediate radial scale-length of 2.5 kpc, and a constant DM density, they then found ρDM = 14 ± 1 mM pc-3. Bienaymé et al. did split their kinematic dataset into three abundance groups, corresponding to [Fe/H] <0.58, [Fe/H] >0.25, and an intermediate population; while not enough to isolate homogenous tracers, this should result in similar vertical scale-heights (see Fig. 4 in Bovy et al. 2012b), so that the density-weighted σz should not be a strong function of height even if each subsample is kinematically inhomogenous. However, interstellar extinction makes the radial range over which these samples were obtained a strong function of galactic latitude: the thin-disc data was obtained over a significantly smaller range of Galactic radii than the thick-disc data. Chen et al. (2012) showed that the thin- and thick-disc scale-lengths estimated from SEGUE stars near the Galactic plane are different, confirming the results of Bensby et al. (2011) using many fewer stars: their values are Hthin = 3.4 and Hthick = 1.8 kpc, suggesting that the thin-disc scale-length is nearly twice that of the thick disc. Assuming that the squared velocity dispersions vary exponentially with the same scale-length and hence are larger at smaller Galactic radii, the thick-disc data is weighted to smaller average radii and hence larger dispersions. Thus, there must be an inverse correlation between the radial scale-length and σz, producing an artificially larger dispersion and hence larger ρDM for the data thought to be the most sensitive. Again, a crude model can be used to estimate the effects: for an object selection criterion roughly corresponding to that used by Bienaymé et al., consisting of a minimum distance of 200 pc, a maximum distance of 3000 pc, | b | > 22°, and radial scale-lengths of 1.6 or 3.4 kpc for the thick and thin disc, respectively, one expects an artificial increase in the relative vertically averaged values of the thick-disc σz over those of the thin disc over the range of 200 to 2000 pc of ~20% independent of the common vertical scale-height. This is roughly the gradient seen in the RAVE data and a sampling effect not considered in their model.

Even when considerable effort is put into an attempt to correct for systematic effects, to treat consistently inhomogenous tracers, and to minimized the number of assumptions, the result is not necessarily more reliable. For instance, Garbari et al. (2011) use a detailed global MW model in an attempt to correct for the effects of a realistic Galaxy (e.g. spiral arms), but their mock-data simulation has neither a thick-disc component nor an ISM component even though these two components make up 6070% of the mid-plane density (Table 1) and differences in radial scale-lengths between the thin and thick discs must also affect the sampling. In addition, their fits were constrained to have ρb,local(0) in the unnecessarily narrow range 91 ± 14 mM pc-3, which implies a considerable shift in kinematic influence to the otherwise less well-constrained DM halo.

The problems with mixed tracer populations can be minimized by carefully choosing each sample to be as homogenous as possible, and the most obvious criterion of homogeneity is detailed abundance, e.g. [Fe/H] and/or [α/Fe]. Zhang et al. (2013; hereafter Z13) derived a DM density of 6.5 ± 2.3 mM pc-3 using three different [α/Fe]-selected, K-dwarf tracer populations taken from the SEGUE survey at vertical heights between 300 and 1400 pc, roughly corresponding to thin-, intermediate-, and thick-disc tracer populations. However, there are several problems with this analysis as well: Z13 used an approximate vertical Jeans equation; they assumed that the baryonic contributions to the local vertical gravitational field are due to infinite homogenous discs; their stellar disc had a thin-disc scale-height even though at the higher heights the gravity should be dominated by the thick disc; they assumed the traditionally low gas surface density (see discussion in the previous section); the three (assumed exponential) tracer populations each cover a wide range in [Fe/H] and hence in scale-heights, in all probability producing the kinematic inconsistencies implied by Eq. (1) (the thin and intermediate tracers have density profiles with non-exponential tails, as expected); and they used an undocumented Markov chain Monte Carlo (MCMC) prior on the scale-heights (Zhang, priv. comm.) without which their value of ρDM is undetermined.

Bovy & Rix (2013) analysed a SEGUE/SDSS dataset of G dwarfs at Galactic radii 4.5 <R< 7 kpc and much larger heights than was possible with the K dwarfs used by Z13. They modelled the data using DF models for each mono-abundance ([Fe/H] and [α/Fe]) data subset and modified Stäckel potentials, permitting them to simultaneously estimate important quantities like the radial scale-lengths. Many publications now use their derived value of the stellar surface density as an assumed accurate measure of the total baryonic contribution (e.g. Iocco et al. 2015). However, the use of Stäckel potentials results in vertical dispersion profiles σz(z) that increase towards the disc mid-plane for heights below 2 kpc, and are otherwise flat (their Fig. 4d); measurements of the near-mid-plane dispersions are significantly lower (Kuijken & Gilmore 1989b; Crézé et al. 1998; Dehnen & Binney 1998; Holmberg et al. 2007; Pasetto et al. 2012a) and the mono-abundance data of Büdenbender et al. (2014; hereafter B14) clearly show that the gradient at such heights, if present, has the opposite sign. Part of this problem may be due to the fact that their model did not account for the effects of the local inter-arm region (see the discussion in the previous section), which would produce an apparent change in the radial scale-length and hence potentially observable changes in the modelled kinematics. They also used the same traditionally low estimate for the gaseous surface density and a single-component disc potential. Other issues are discussed by Read (2014). Thus, while Bovy & Rix (2013) have shown how potentially powerful the complicated method for modelling the kinematics of disc stars can be, their impressively precise results are, in fact, inaccurate at a level that is important when trying to measure ρDM.

Büdenbender et al. (2014) separated the SEGUE G dwarf data into clearly thin- and thick-disc tracers and showed that the simple model used by K13 is unable to simultaneously explain the differences in the magnitudes and the slopes of the two tracers. They attributed this to different orientations of the velocity ellipsoid with height, implying that the vertical motions are not fully decoupled from the radial ones. However, both the thin- and thick-disc tracer populations were not mono-abundance selected and, like the Z13 samples, contained stars with very different [Fe/H] and hence vertical scale-heights. Was the difference due to model assumptions and fundamental kinematic effects, or simply due to inhomogeneous tracers?

thumbnail Fig. 1

The vertical gravity gz of axisymmetric doubly exponential MW discs and radial scale-lengths H = 1.4,1.8,2.2,2.6,3.0,3.4, and 3.8 kpc; the blue line is gz for a uniform infinite plane with the default local surface density Σ0 ≡ Σ(R0;A ≡ 0); the red lines are fits to gz(z< 3 kpc) assuming an additional effective infinite uniform disc component (see text).

Open with DEXTER

4. The not-so-local gravity field

In order to derive ρDM from homogenous velocity dispersion data using a minimally complicated model of the MW – stationary, axisymmetric, vertically symmetric – in a manner similar to that used by Z13, one must specify the gravitational forces on the stars more accurately than possible in the simple model they used. Although the difference between the gravity from an infinite, uniform disc and the real MW should be small for | z | ≪ R0, the MBCM data reach to ~R0/ 2 and a correction to the model of an infinite uniform disc that takes structure at larger scales is needed. Unfortunately, the effort to construct a global mass model is exactly what one wanted to avoid by solving the vertical Jeans equation alone, so the question is whether the non-uniform effects are significant and, if so, how easily they can be included as corrections to a local analysis.

The effects due to the global variation of the disc surface density with radius, parameterised by the disc scale-length H of some component, are straightforward to calculate and are often – but by no means always – included in kinematical analyses when they should be. For large H, the disc regions at smaller and larger Galactic radii will have similar surface densities and one would naively expect the correction to be small. For small H, the surface density at smaller Galactic radii will have a significantly larger surface density that could dominate the local vertical gravity contribution and which may or may not be compensated by the loss of surface density at larger radii – the gravitational field of a planar annulus is very different from that of a point source, making it hard to know the outcome a priori. For a doubly exponential disc component, this effect can be calculated semi-analytically using the expressions for gz(R,z) derived by Kuijken & Gilmore (1989c)3.

The results for a doubly exponential disc with vertical scale-height h = 300 pc and a range of H are shown in Fig. 1 (top): for H ≈ 3 kpc, gz is very similar to that due to an infinite uniform disc; for larger H there is a slight decrease at large | z | ≫ h; and for smaller H there is a substantial increase in magnitude and the lack of a flattening at large | z | ≫ h. Note that the last effects increase the contribution to gz from the baryonic discs – perhaps dramatically – compared to a model that uses local disc properties only, and the second one mimics the effects of a constant DM density (gz(disc) ≉ const. as | z | → ∞); both must result in a reduction in the amount of DM needed to fit the data. This is the opposite effect to that reported by Kuijken & Gilmore (1989b): in their Fig. 2, they assumed a smaller R0 = 7.4 kpc and a larger H = 4.5 kpc, i.e. H/R0 = 0.6. For values of H/R0 near unity, the gravitational field of the disc starts to approach that of a point mass, leading to the drop in gz(z) at large heights.

The change in the form of gz relative to that of an infinite disc is itself very similar to that of an infinite uniform disc with a different surface density and scale-height: in Fig. 1, the fits to the real form of gz(z) consists of that for an infinite disc with the nominal local surface density plus an additional infinite disc component with fitted surface density and scale-height (e.g. the decrease in gz(z) at large | z | for large H is then modelled as a negative additional surface density). When the Jurić et al. (2008) thin- and thick-disc vertical scale-heights h of 300 and 900 pc and the corresponding radial scale-lengths H from Bensby et al. (2011) of 3.4 and 1.8 kpc, respectively, are used, the gravity due to the thin disc is nearly unaffected whereas the effective thick-disc surface density increases by 140%, resulting in an increase in the total effective surface density of 37%, and the additional component has a large effective vertical scale-length of 2.6 kpc. The changes in the form of gz(z) for the very thin gaseous disc of the MW is negligible because the radial scale-length is roughly twice that of the disc (Bigiel & Blitz 2012).

Thus, the effects on the vertical gravity profile of non-uniform stellar structures in the Galactic neighbourhood of the Sun are quite significant and must be included in all vertical Jeans analyses. The naive separation of the gravity into an asymptotically constant disc and increasing DM component (e.g. the classic Kuijken & Gilmore “K” and “F” models) is not necessarily an adequate representation. Fortunately, even though the local vertical gravity is not that of an infinite disc with constant surface density, the effects of a realistic axisymmetric mass-distribution for the disc of the MW can be modelled using additional infinite disc components whose surface densities and scale-heights depend simply upon the true values.

thumbnail Fig. 2

and data from Büdenbender et al. (2014; B14) and Pasetto et al. (2012a,b); for the B14 SEGUE data, each colour represents a different [Fe/H] and [α/Fe] mono-abundance tracer population (and was chosen to reproduce the corresponding plots in B14), whereas the Pasetto et al. data are for undifferentiated RAVE stars.

Open with DEXTER

5. The vertical Jeans equation revisited

The vertical Jeans equation can be easily derived from moments of the collisionless Boltzmann equation in cylindrical coordinates for a steady, axisymmetric stellar tracer component i with density ρi(R,z) (see BT08; Eq. (4.222b), p. 353). For a steady-state situation ( for k = R,z, i.e. no net flux of stars in R or z), one obtains the vertical, axisymmetric Jeans equation (2)This is an exact and full description of the connection between the local vertical gravity and these local stellar velocity dispersions within the very modest assumptions made. Admittedly, there is evidence that the vertical structure of the MW’s thin disc is not stationary (e.g. Widrow et al. 2012; Carlin et al. 2014), but these are second-order effects that should not seriously affect the velocity dispersions, particularly at larger heights. Admittedly, the gravitational field of the MW is not axisymmetric owing to factors like the central bar and the spiral arms, but the effects of the former are subtle (e.g. Bovy et al. 2014) and the non-axisymmetric effects of the latter are small due to local pitch angles of less than 20° (Vallée 2005; Levine et al. 2006).

Fortunately, the mean radial surface brightnesses of disc stars in most spiral galaxies are known to be exponential (Freeman 1970) and the local vertical distributions of stars in the local MW are also observed to be exponential starting at quite modest heights above the plane of the disc (Jurić et al. 2008). In particular, this holds for mono-abundance tracers (Bovy et al. 2012a,b; K13) so that one can confidently assume a doubly exponential density profile for any homogenous tracer population, and simply measure the radial scale-lengths Hi and vertical scale-heights hi in order to characterise this most important property to an adequate accuracy : (3)Given the observed properties of ρi, it is not unreasonable to make a similar assumption about the radial variation of the velocity dispersion: (4)A simple exponential disc with a constant disc thickness h requires σzz ∝ Σ, i.e. LzH. More complicated dynamical models of the MW suggest that Lz is larger than H by ~20% (e.g. Sanders & Binney 2015), so assuming they are equal should not be a major source of error.

To solve the vertical Jeans equation fully, one needs to know the radial and vertical behaviour of σRz. This component generally has a smaller effect than σzz because the hi are generally much smaller than the Hi, so this contribution was totally left out in Z13. Several possible models for σRz have been proposed (see Amendt & Cuddeford 1991). For example, one can assume that the principle axes of the velocity ellipsoid are aligned in spherical or cylindrical coordinates so that the correlation between the radial and vertical motions can be expressed via a tilt angle αtilt such that (5)(see BT08, problem 4.32, p. 391). For the special case of Stäckel potentials, it can be shown that αtilt is aligned with the potential gradients and so should be independent of the tracer population used (e.g. Statler 1989; Binney 2012). This form was used by Siebert et al. (2008) and Cassetti-Dinescu et al. (2011) to analyse abundance-undifferentiated data from the (southern) RAVE experiment close to the Galactic plane. Büdenbender et al. (2014) have measured αtilt for the abundance-differentiated dwarf G stars in the (northern) SEGUE sample and saw that αtilt increases in magnitude with increasing z: they fitted αtilt as a linear function in z but, given the errors, an equally good fit to the data is simply (6)Both fits imply that the tilt angle is essentially aligned with the local spherical coordinate system. If an anisotropy parameter is defined (e.g. Pasetto et al. 2012b), the assumptions β ≡ const. plus a tilt angle perfectly aligned with the Galactic center result in (7)(Kuijken & Gilmore 1989c), implying that the assumption c ≈ const. is not unreasonable for z2βR2.

The SEGUE data from the B14 Table 1 are shown in Fig. 2, covering heights 0.5 <z< 2.5 kpc. Clearly, there is a tight connection between and that appears to be quite linear within a given tracer population and very similar for different tracer populations4. If one generically assumes (8)(i.e. βi ≡ 1 /bi for ai ≡ 0), one obtains a single equation expressing σRz in terms of σz, z, and two dispersion-model parameters η and γ, (9)where η ≡ −ca/ 2 and γ ≡ −c(b−1)/2 (the sign was chosen to fit the northern SEGUE data of B14). While this expression was derived using the tilt angle formalism, it constitutes a simple, phenomenological model for σRz in terms of the model parameters γ and η and an assumed scaling in z that could cover other origins for the behaviour of this dispersion component. The relatively simple assumption that σRz is proportional to a linear function in σzz whose constant term might be non-zero should not be construed as a potential physical problem: the fact that the function is non-zero for σzz = 0 is irrelevant given that the dispersions are nearly isothermal for any given tracer population, and that what is needed is merely a functional connection between the two dispersions measures.

The usefulness of this model for σRz depends, of course, upon how well it fits the data, at least locally. In order to augment the B14 data with observations at lower heights, the thin-disc RAVE results from Pasetto et al. (2012a) were modelled using Eq. (4) for σzz and the epicyclic approximation for σRR (Amendt & Cuddeford 1991). The fits to σkk(R0,z) were made using the python MCMC object emcee (Foreman-Mackey et al. 2013) with the parameters Lk and the σk values at each height bin (0, ±200, and ± 400 pc). A standard probability function lnP ∝ −χ2 and 30 MCMC walkers were used; after a burn-in of 2000 iterations each, the walkers were well converged and the probability density functions (PDF) were calculated from a total of 150 000 samples. While the values of σzz(R0,z) and σRR(R0,z) are reasonably well constrained, the final PDFs of the scale-lengths have long tails; if the radial scale-lengths are restricted to be with the reasonable range 0 <Lk< 10 kpc, , and (68.2% bounds), consistent with being equal to a thin-disc density scale-length of ~3 kpc. When the RAVE dispersion points are added to Fig. 2, it appears that η ≡ 0, β ≈ const. is a reasonable assumption for the thin disc, but the B14 thick-disc data shows a much flatter relation. Thus, Eq. (9) appears to be an adequate – if not excellent – local phenomenological model that is able to model the potentially different behaviours of different tracers.

With these assumptions, one can then rewrite the local vertical Jeans equation as the ordinary differential equation (10)where (11)is a typical scale-length where the effects of cross-talk between σRz and σzz appear, and λi/hi is a measure of the effects of the σRz term on σzz for a given tracer: larger values imply less effect. The ηi term in Eq. (10) is not proportional to σzz and so behaves like a gravity term. Indeed, the proportionality to z makes it look like a constant apparent mass density defined by (12)which represents a small but tracer-dependent and hence potentially noticeable correction to any fitted DM density for reasonable values of ηi and Hi.

The equation for gz(R0,z) used by Z13 to close the equation assumes an infinite uniform disc plus a constant DM density: in Sect. 4, we saw that a generalisation of their model – the sum of effective infinite discs with different effective surface densities Σj and effective exponential vertical scale-heights lj (the latter chosen to be easily distinguishable from the tracer scale-heights hi), plus a constant DM density ρDM(R0,0) – is a reasonable functional model: (13)This assumption also permits one to solve the vertical Jeans equation (Eq. (10)) analytically when the simple model for σRz (Eqs. (9) and (12)) is used, (14)where

(15)and where (16)is Dawson’s function (erfi(x) is the imaginary error function). This solution is only good for reasonable heights , but this is not a critical constraint on the applicability of the solution since λi/hi ≫ 1 for any reasonable value of γi and hi (Eq. (10)) and the whole procedure breaks down (e.g. the assumptions of exponential density profiles and constant DM densities) if carried out too far from the disc plane.

The behaviour of this solution relative to that of Z13 is governed by the ratio λi/hi, as is shown in Fig. 3. The effects of λi are just apparent for λi/hi ~ 16, but are quite significant for smaller values. Thus, γi can significantly change the form of the solution. In the limit λi → ∞, i.e. γi → 0, this solution is the same as that used by Z13 and B14. However, for small λi/hi, the simple asymptotic behaviour of the Z13 solution (flat disc and linear DM contributions to ) is modified so that both mass components produce increasing σz with height, making both the baryons and the DM more effective at explaining positive dispersion gradients and thus principally reducing the total amount of DM needed.

thumbnail Fig. 3

Top: the function representing the effect of a vertically exponential disc on σz for the ratio of vertical scale-heights lj/hi = 0.5 (red), 1 (blue), and 2 (black) for scaled σRz parameters λi/hi = 4 (highest), 8, 16, and 32 (lowest). Bottom: the same for the function representing the effect of a constant density.

Open with DEXTER

6. Thin-disc versus thick-disc kinematics

Büdenbender et al. (2014) showed that the σz(R0,z) profiles of the thin- and thick-disc tracers in the SEGUE G dwarf dataset look very different (their Fig. 2): the former show a large gradient of ~6 km s-1 kpc-1 at relatively low Galactic heights (| z | < 1.3 kpc) and the latter half the gradient at much larger distances. As emphasised by B14, the simple model for σz used by Z13 is incapable of simultaneously explaining both the relative magnitudes and the shapes of the two curves: one expects a rapid rise in at heights z<hdisc to a constant value due to the disc and a subsequent linear rise due to the constant DM density (Eq. (1)). The difference between different tracers is a simple scaling with the tracer scale-heights hi at large heights hi>lj, so the relative magnitudes are fixed by the relatively well-determined vertical density profiles; they cannot be changed at will, for example by modifying the relative amounts of baryonic and DM. Given that the full vertical Jeans equation can now be used and that one can model the local vertical gravity more realistically, is it now possible to explain these differences, or are they simply due to tracer inhomogeneities?

At first glance, the new model would appear to be able to solve the problem entirely: different values of γi (or, equivalently, of λi) permit either a flat or an increasing asymptotic σz profile in the same disc but for different tracer populations, exactly as seen by B14. The relative magnitudes are affected by different values of either γi or ηi, and by the difference due to relative values of the vertical scale-heights hi. In contrast to the Z13 model, a flat profile does not necessarily require unrealistic amounts of baryons, since a small exponential radial scale-length can increase the effective surface density used to calculate the gravity by large amounts. If we find the right values for the kinematic and gravity model parameters, there would appear to be nothing in the way of determining a robust value for the DM density.

As we saw in Sect. 4, a short radial scale-length H ≈ 2 kpc for the thick-disc means that the importance and effective scale-height of this gravity component is much larger than previously considered, placing its effects in a region where an increase in σz was thought to be due to DM alone (Fig. 1). Since the gravitational effects of the thick-disc should gradually increase with height, whereas those of the thin disc eventually saturate, there should be an additional gradient in σz for the thick-disc tracer data which samples the upper regions. This is, unfortunately, exactly the opposite of what is seen in B14’s two tracer samples, so the explanation is not simply a question of the gravity model.

The main kinematic parameters are ηi and γi, expressible as the effective mass-density (Eq. (12)) and the scale-length λi (Eq. (14)). The value of is proportional to Hi/γi, so the smaller the value, the more significant the effects. While B14 did not publish σRz or σRR values for their two representative thin- and thick-disc samples, one can take the kinematic data in their Table 1 and split them into two kinematic groups, corresponding to the clearly thin discs ([α/Fe] < 0.2) and clearly thick discs ([α/Fe] > 0.28), i.e. the intermediate disc tracer with [Fe/H] = 0.35 and [α/Fe] = 0.28 is excluded. This makes it possible to attempt to fit σRz and hence attempt to explain the differences seen by B14.

The first model (“A”) uses the tilt-angle formalism connecting σRz with both σR and σz (Eq. (5)). This entails fitting the model parameters a and b describing σR(σz) (Eq. (8)) and then using these results to fit σRz(z,σR,σz) via the parameter c, from which γ, η, and λ can finally be derived for the B14 thin, thick, and combined datasets (model “A1”). The fits were again made using emcee (Foreman-Mackey et al. 2013) following a similar fitting process to that described in Sect. 2. The value of for the thin dataset is large (3.1 vs. 1.1 for the thick dataset), but the parameter c representing the tilt-angle behaviour is clearly different for the two populations (−1.4 ± 0.2 vs. −2.5 ± 0.3) unlike what one would expect for a gravitational potential representable by Stäckel functions. The ratio λ/h is also different (13 ± 3 vs. 3.5 ± 0.5), suggesting that the thick-disc σz should have a larger gradient than the thin disc, again exactly the opposite of what is seen in B14. Assuming that the tilt-angle behaviour is common to both tracers, one can fit a common c parameter (model “A2”); this fit shows similar values for the model parameters and the same unexpected λ/h behaviour between the thin- and thick-disc tracers. Since there is no reason why the σR(σz)-relation should be the same for all tracers, one can use an extreme model which avoids the use of η by setting a and η equal to zero; this model (“A3”) provides adequate explanations for σR(z,σz) when applied to the thin- and thick-disc data separately but results in very bad fits for the combined data (). Despite this drastic modification of the kinematic model, the same unexpected λ/h behaviour is seen.

The alternative phenomenological model (“B1”) is to adopt the mathematical form of σRz(z,σz), but to fit the two model parameters γ and η directly, ignoring the question of the implied effect on σR(σz) within the tilt-angle assumption. The values of for this model are lower simply because we are now fitting only the σRz data rather than the higher quality σR data as well, but the fitted values of η are not significant and the same λ/h behaviour is seen. Finally, one can fit a one-parameter model (“B2”) by assuming η ≡ 0; these fits are just as good ( is 1.2 or better) and σRz is well fit, but one sees the same λ/h behaviour: γ is 1.3 ± 0.2 and 0.8 ± 0.1, and λ/h is 8.9 ± 0.8 and 3.2 ± 0.2 for the thin- and thick-discs, respectively.

Thus, one must conclude that the differences in the kinematic behaviours of the B14 thin- and thick-disc samples are not simply due to the effects of the cross-dispersion term σRz, at least as modelled using the tilt-angle paradigm. The fact that a realistic gravity model should also result in the opposite behaviour seen by B14 reinforces this conclusion. If the effects seen in the SEGUE thin- and thick-disc data by B14 are not simply due to the inadequacies of the kinematic model, they must be due to the data. For example, the larger gradient in the B14 σz thin-disc data is mainly due to the outermost z point at 1.25 kpc: if this point is lowered, then the rise at lower heights can be due to the disc scale-heights alone. Alternatively, the effect might be due to the same behaviour seen in the MBCM data – a gradient produced by mixing the different kinematic behaviours of different tracers. If one compares the B14 sampling with that of Bovy et al. (2012a, their Fig. 4), one sees that the vertical scale-heights of the thin-disc mono-abundance tracers in the B14 thin-disc samples vary from about 200 to 400 kpc: using the same simple calculation made for the BT12 analysis of the MBCM data in Sect. 3, one can produce an increase from 18 to 25 km s-1 over the height range of 0.40 to 1.25 kpc (Fig. 3 in B14) with a local density ratio between the thinest and the least thin isothermal tracer component because the ratio of the relative contributions of the former over this range is roughly 10% at the upper regions of the thin-disc data. For a reasonable mixture of densities, one then expects an artificial gradient even in the thin-disc σz data due solely to population mixture. A similarly strong effect for the thick-disc data is not expected, since the thick-disc sample is dominated by the two most extreme thick-disc tracer populations which show roughly the same asymptotic values of σz, but an artificial increase in the inhomogenous vertical dispersion data must be present here as well.

If there is no problem in the vertical gradients of σz measured for truly homogenous data, is there one in the relative magnitudes, as suggested by B14 (the dashed solution in their left Fig. 2 diagram)? The B14 “no dark matter” solution used a very large baryonic surface density to explain the thick-disc data, but the model overpredicted the σz level for the thin disc by about 20% because the asymptotic ratio is simply hthick/hthin for the Z13 model, independent of the relative contributions of baryonic discs and the DM halo. In the model for σz presented here, the ratio also depends upon the relative values of λi/hi and ηi, i.e. in a complicated fashion that depends upon the full details of the models and the data used. A full discussion of these possibilities must be relegated to a subsequent paper, but it suffices to say that critical uncertainties in the vertical scale-heights and the obviously different values of γi are enough to reduce this discrepancy to a level where the model is principally able to explain the observations. Thus, the simplest explanation for the thin- and thick-disc kinematic behaviour seen by B14 is that the tracer data binned to represent only a thin and a thick disc is not homogeneous enough to permit the simple kinematic analysis used.

7. Conclusions

I have shown that many previous attempts at estimating the local DM density have not been as robust as traditionally perceived. The usual baryonic corrections to the local dynamic density underestimate or ignore the uncertainties in the traditional gaseous densities and the traditional densities did not include the substantial amounts of “dark” gas in the form of optically thick HI and CO-dark molecular gas now known to be present in the ISM. The non-local effects of Galactic structure increase the effective dynamic density by another 2535%, resulting in a total dynamic baryonic density of 120−160 mM pc-3, leaving practically no room for measuring DM using mid-plane kinematics alone. The mixing of different tracer populations with different vertical scale-heights, radial scale-lengths, and dispersion magnitudes produces gradients in σz(z) that can be misinterpreted as being due to DM whenever the effects of kinematic inhomogeneity are not properly taken into account. The gravitational effects of doubly exponential discs are quite different in form and magnitude from the infinite discs often used, particularly for radial scale-lengths H< 3 kpc. Fortunately, one can easily model this effect using additional infinite disc components with fitted surface densities and scale-heights. For a simple model of the MW thin and thick discs, the increase in the effective total gravitating surface density is quite large, ~40%.

Assuming only axisymmetry, stationarity, and planar symmetry, it is possible to derive a full analytic solution of the vertical Jeans equation for homogeneous tracers, a reasonably realistic local MW disc plus a constant DM density that removes many of the approximations present in previous analyses. The dispersion cross-term σRz can be been modelled using a phenomenologically justified expression motivated by, but not dependent upon, a tilt-angle assumption. The additional dependence on σz in the vertical Jeans equation can induce tracer-dependent and detectable additional vertical dispersion gradients that can systematically reduce the amount of both baryonic and dark matter needed to explain the observations.

The σRz model parameters were estimated using the SEQUE G dwarf data analysed by B12 and B14, either by assuming a tilt-angle model or by fitting a simple phenomenological model dependent on a linear function in and a factor of z/R. Despite the additional complexity of form available, it is not possible to explain the different behaviour of the “thin” and “thick” tracer data of B14 using kinematical effects alone: while the cross-dispersion effects should be large enough to be detected in the vertical profile of σz(R0,z), the changes are opposite to those needed to explain the different dispersion behaviours in the two B14 samples. The simplest explanation is that these two inhomogeneous (non-mono-abundance) samples also show the tracer population mixture effects plaguing other analyses.

Thus, an analysis of local stellar dispersion data to extract a robust local value of ρDM should in fact be possible, but it would require considerably more care and better data than has been used to date (e.g. gaia). In a subsequent paper, I will apply the techniques outlined herein to carefully chosen data in order to estimate the local DM density.


1

I use the astronomical unit 1 mM pc-3 ≡ 0.001 M pc-3 = 0.038 GeV cm-3 to avoid having to count the many zeros present when using M pc-3.

2

Read (2014) forgot to correct the Kalberla & Dedes value of ~12 M pc-2 from the pure hydrogen to the total density.

3

The vertical acceleration gz is often labelled Kz; the latter often has an ill-defined sign, and is often incorrectly called a “force” in the literature.

4

Actually, all of the data are consistent with the single assumption , but there is no obvious reason for this relation and it is not very useful for the analysis that follows.

Acknowledgments

I would like to thank J. Bovy, A. Büdenbender, and L. Zhang for kind access to their data. I particularly appreciate detailed comments by C. Moni Bidin and the anonymous referee for their helpful comments that led to considerable improvements in this manuscript. This work used the HI galaxy images available at http://www.mpia-hd.mpg.de/THINGS/Overview.html courtesy of the THINGS collaboration.

References

All Tables

Table 1

A local baryonic mass model using classic (updated) ISM density estimates.

All Figures

thumbnail Fig. 1

The vertical gravity gz of axisymmetric doubly exponential MW discs and radial scale-lengths H = 1.4,1.8,2.2,2.6,3.0,3.4, and 3.8 kpc; the blue line is gz for a uniform infinite plane with the default local surface density Σ0 ≡ Σ(R0;A ≡ 0); the red lines are fits to gz(z< 3 kpc) assuming an additional effective infinite uniform disc component (see text).

Open with DEXTER
In the text
thumbnail Fig. 2

and data from Büdenbender et al. (2014; B14) and Pasetto et al. (2012a,b); for the B14 SEGUE data, each colour represents a different [Fe/H] and [α/Fe] mono-abundance tracer population (and was chosen to reproduce the corresponding plots in B14), whereas the Pasetto et al. data are for undifferentiated RAVE stars.

Open with DEXTER
In the text
thumbnail Fig. 3

Top: the function representing the effect of a vertically exponential disc on σz for the ratio of vertical scale-heights lj/hi = 0.5 (red), 1 (blue), and 2 (black) for scaled σRz parameters λi/hi = 4 (highest), 8, 16, and 32 (lowest). Bottom: the same for the function representing the effect of a constant density.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.