A non-parametric method for measuring the local dark matter density

Silverwood, H.; Sivertsson, S.; Steger, P.; Read, J. I.; Bertone, G.

doi:10.1093/mnras/stw917

Abstract

We present a new method for determining the local dark matter density using kinematic data for a population of tracer stars. The Jeans equation in the z-direction is integrated to yield an equation that gives the velocity dispersion as a function of the total mass density, tracer density, and the ‘tilt’ term that describes the coupling of vertical and radial motions. We then fit a dark matter mass profile to tracer density and velocity dispersion data to derive credible regions on the vertical dark matter density profile. Our method avoids numerical differentiation, leading to lower numerical noise, and is able to deal with the tilt term while remaining one dimensional. In this study we present the method and perform initial tests on idealized mock data. We also demonstrate the importance of dealing with the tilt term for tracers that sample ≳1 kpc above the disc plane. If ignored, this results in a systematic underestimation of the dark matter density.

Galaxy: disc, Galaxy: kinematics and dynamics, dark matter

1 INTRODUCTION

Dark matter (DM) is an elusive component of the cosmos. Its existence has long been inferred from its gravitational interactions with ordinary matter on scales ranging from dwarf galaxies and galactic rotation curves to lensing galaxies, galaxy clusters, and the Universe as a whole (for reviews see e.g. Rubin 1993; Jungman, Kamionkowski & Griest 1996; Bergström 2000; Bertone, Hooper & Silk 2005; Bertone 2010; Massey, Kitching & Richard 2010). However its exact nature remains unknown.

Understanding the distribution of DM in cosmological structures, and in particular in the Milky Way, is of great importance to test the standard Λcold dark matter (ΛCDM) cosmological model, to make predictions – and allow the interpretation – of experiments that seek to detect DM, and to probe more exotic models for DM (e.g. Read 2014). Here, we will focus our attention on the local DM density – i.e. the average density of DM over a small volume (typically of size ∼100–1000 pc) centred around the Sun. The local DM density is a vital input for so-called direct DM experiments, based on the detection of the recoil energy of nuclei struck by DM particles streaming through underground detectors, as the expected event rate is proportional to the local DM density and the DM–nucleon cross-section. In addition, the local DM density is important for indirect searches that look for neutrinos produced by the annihilation of DM particles trapped at the centre of Sun (see e.g. Klasen, Pohl & Sigl 2015, and references therein for a recent update on direct and indirect searches). The results of these searches are also used in explorations of theoretical parameter space, such as those of Supersymmetry (e.g. Allanach & Lester 2006; Buchmueller et al. 2014; Strege et al. 2014; Bertone et al. 2016; Roszkowski, Sessolo & Williams 2015).

Methods of determining the local DM density can be divided into two categories: those utilizing measurements of stars in a small volume around the Sun (e.g. Kapteyn 1922; Jeans 1922; Oort 1932; Bahcall 1984; Kuijken & Gilmore 1989b, 1991; Creze et al. 1998; Garbari et al. 2012; Bovy & Tremaine 2012; Smith, Whiteoak & Evans 2012; Zhang et al. 2013; Bienaymé et al. 2014) and those utilizing a host of dynamical tracers – particularly the rotation curve – to constrain more global mass models of the Milky Way (Dehnen & Binney 1998; Catena & Ullio 2010; Salucci et al. 2010; Weber & de Boer 2010; McMillan 2011; Nesti & Salucci 2013; Piffl et al. 2014; Iocco, Pato & Bertone 2015; Pato & Iocco 2015; Pato, Iocco & Bertone 2015). Following the notation of Read (2014) we denote the local DM density derived from local measurements by ρ_DM, and those extrapolated from rotation curves assuming a spherical DM halo as ρ_{DM, ext}. With the advent of large survey data, local and global methods are beginning to converge (Bovy & Rix 2013; Piffl et al. 2014).

The comparison of ρ_DM to ρ_{DM, ext} can provide insight into the shape of the Milky Way's DM halo, and thus into the formation of the galaxy (Read 2014). If ρ_DM from local measurements is smaller than that extrapolated from global measurements, ρ_{DM, ext}, then this implies a prolate halo, while the opposite would imply a squashed, oblate halo. While the former is produced by DM only N-body simulations, the introduction of baryons into the simulations produces the latter (Katz & Gunn 1991; Dubinski 1994; Debattista et al. 2008; Read et al. 2009).

If ρ_DM is derived from tracers distributed vertically in the plane it is possible to probe not only the local DM density but also its vertical profile, ρ_DM(z). This could potentially reveal deviations from a spherical halo profile, which can form in a number of ways. First, the DM halo should respond to the formation of the baryonic disc by flattening into an oblate halo (Katz & Gunn 1991; Dubinski 1994; Debattista et al. 2008; Read et al. 2009). Secondly, the accretion of subhaloes after the formation of the baryonic disc can also lead to the formation of a ‘dark disc’ (DD) that corotates with the baryonic disc (Lake 1989; Read et al. 2008, 2009; Purcell, Bullock & Kaplinghat 2009; Ling et al. 2010; Pillepich et al. 2014). Gravitationally these two features, an oblate halo and a DD, are indistinguishable and degenerate. One method to distinguish the two is to hunt for the chemically and kinematically distinct sub-halo stars that would accompany the accreted DD, but not the contracted DM halo. Recent Gaia-ESO observational data find no evidence for such stars (Ruchti et al. 2014, 2015), suggesting that the Milky Way has had a rather quiet past since its disc formed, and therefore does not host a significant accreted DD. If correct, then any gravitationally detected non-sphericity must imply a non-rotating locally oblate halo. An enhancement in ρ_DM over ρ_{DM, ext} is also a feature in some more exotic models, such as Partially Interacting Dark Matter (PIDM) (Fan et al. 2013), and some modified theories of gravity (e.g. Read & Moore 2005; Nipoti et al. 2007).

Measurement of the local DM density dates back almost 100 years (Kapteyn 1922). The history of this measurement is one of both increasingly precise data and decreasingly strong assumptions. As such, the error bars on ρ_DM have not always shrunk with time, as better data often allows one to discard previous assumptions. With the advent of Gaia data from 2016 onwards we will have access to high-precision data on individual stars, and the primary uncertainty in the determination of ρ_DM will be systematic model uncertainties. Recent results for the local DM density include Bienaymé et al. (2014), who used RAVE data to derive a value of ρ_DM = 0.0143 ± 0.0011 M_⊙ pc⁻³ = 0.542 ± 0.042 GeVcm⁻³, and Zhang et al. (2013) who used Sloan Digital Sky Survey (SDSS)/SEGUE data to find ρ_DM = 0.0065 ± 0.0023 M_⊙ pc⁻³ = 0.25 ± 0.09 GeVcm⁻³. Note that these two results do not overlap within their stated uncertainties. The significance of this discrepancy is hard to interpret though, as they each use different data sets and analysis techniques, and both methodologies rely on rather different assumptions. To make progress we should endeavour to minimize the number of assumptions made, and apply the same analysis techniques to both data sets.

In this paper we make progress towards the first goal of reducing model assumptions. To achieve this we introduce a one-dimensional Jeans analysis method to probe the local DM density using the vertical motions of tracer stars. We construct a representation of the tracer density ν, and also allow for a DM density profile that is more complex than simply constant with height as previous local studies have assumed. Additionally, we deal with the so-called tilt term, which links radial and vertical motions of the tracer stars. This term is crucial to understand stellar motions at high-z, where the baryonic contribution falls off and DM becomes increasingly important. Using the vertical Jeans equation we calculate the velocity dispersion σ_z(z) for each mass model, and then fit to data in ν, σ_z, and σ_Rz using MultiNest (Feroz & Hobson 2008; Feroz, Hobson & Bridges 2009; Feroz et al. 2013). We test this method on mock data sets, and explore the impact of tracer star sample size, the tilt term, and non-constant DM profiles. The mock data for this paper are ‘as good as it gets’ – we assume the population to have no observational biases and there to be no measurement error on individual stars – allowing us to isolate the effects of sampling error and model uncertainties. This is similar to what we will have with Gaia data, where the measurement errors on stars will be small compared to sampling error. We will explore the effect of adding realistic Gaia uncertainties to our method in a separate work.

We reiterate that this method is one dimensional, allowing us to keep assumptions to a minimum. However, we show that we are still able to deal with tilt and high-z data, which usually require a two-dimensional method. The key disadvantages of our method are that we bin data and thus lose information, and that our strictly local method cannot take advantage of the ‘long lever arm’ of a global model that would ensure, for example, that the DM density in radial slices close to the ‘local volume’ is continuous and smooth. However, the effect of these disadvantages is to overestimate the errors on ρ_DM, which is acceptable as we aim for a conservative and robust estimate of ρ_DM.

The paper is organized as follows: in Section 2 we introduce our method, covering the Jeans equation based mathematical formalism, our treatment of the rotation curve and tilt term, our descriptions of the elements of our mass and tracer density models. We then outline our statistical analysis, introducing the framework for Bayesian parameter estimation and the MultiNest nested sampling code. In Section 3 we describe the array of mock data sets we use to test out methods. In Section 4 we present and discuss the results of our analyses, before finally concluding the paper with Section 5.

2 METHOD

The broad picture of our problem is this: we have quantities derived from the motions of the tracer stars, namely ν, the tracer density, σ_z, the vertical velocity dispersion, and σ_{R, z}, the (R, z) velocity dispersion. Then we have elements of the mass profile, one of which is unknown – the DM density ρ_DM, and the other which is known within a band of uncertainty – the baryon density ρ_baryon. In this section we first derive the equations to link these quantities, and then describe how each is modelled.

2.1 Deriving a general 1D jeans method

A population of ‘tracer’ stars moving vertically near the Sun obeys the collisionless Boltzmann equation:

\begin{equation} \frac{\mathrm{d}f}{\mathrm{d}t} = \frac{\mathrm{\partial} f}{\mathrm{\partial} t} + \nabla _x f\cdot {\boldsymbol v} - \nabla _v f\cdot \nabla _x \Phi = 0 \end{equation}

(1)

where |$f({\boldsymbol x},{\boldsymbol v})$| is the stellar distribution function; |${\boldsymbol x}$| and |${\boldsymbol v}$| are the positions and velocities, respectively; and Φ is the gravitational potential.

If the system is in dynamic equilibrium (steady state), then we may neglect the partial time derivative of f. Thus for equilibrium tracers, if we know their phase space distribution function f, then we can directly measure the gravitational force from equation (1). In practice, however, this is hard because f is six-dimensional (even a million stars give only 10 sample points per dimension) and we need to estimate the partial derivatives of f to solve equation (1). For this reason, it is common to integrate equation (1) over velocity to obtain a set of moment equations: the Jeans equations (Binney & Tremaine 2008). Adopting cylindrical polar coordinates (R, ϕ, z) and focusing on the z-Jeans equation, we have

\begin{equation} \underbrace{\frac{1}{R\nu }\frac{\mathrm{\partial} }{\mathrm{\partial} R}\left(R\nu \sigma _{Rz}\right)}_{{\rm tilt \,term}:\,\,\mathcal {T}} + \underbrace{\frac{1}{R\nu }\frac{\mathrm{\partial} }{\mathrm{\partial} \phi }\left(\nu \sigma _{\phi z}\right)}_{{\rm axial {'} \, term}:\,\,\mathcal {A}} + \frac{1}{\nu }\frac{\mathrm{d}}{\mathrm{d}z}\left(\nu \sigma _z^2\right) = \underbrace{-\frac{\mathrm{d}\Phi }{\mathrm{d}z}}_{K_z} \end{equation}

(2)

where ν and |$\sigma _z^2$| are the density and vertical velocity dispersion profile of a tracer population as a function of height z above the disc plane, σ_Rz is the cross term in the velocity dispersion tensor that couples radial and vertical motions, and σ_ϕz is the cross-term coupling vertical and axial motions.

Integrating both sides of equation (2), we derive our key equation for this paper:

\begin{equation} \sigma _z^2(z) = \frac{1}{\nu (z)}\int _{0}^z \nu (z^{\prime }) \left[K_z(z^{\prime }) - \mathcal {T}(z^{\prime }) - \mathcal {A}(z^{\prime }) \right] {\rm d}z^{\prime } + \frac{C}{\nu (z)} \end{equation}

(3)

where C is a normalization parameter. For z = 0, we have that

\begin{equation} \sigma _z^2(0)\nu (0) = C, \end{equation}

(4)

and so C simply sets the vertical velocity dispersion at z = 0. [Note that a similar but less general equation was derived recently in Smith et al. (2012), equation (10).] We implement the constant C as a parameter in our model that we will ultimately marginalize over. The alternative would be to calculate C directly from equation (3), utilizing the fact that as z → ∞, ν(z) → 0, which gives us

\begin{equation} C =-\int _{0}^\infty \nu (z^{\prime }) \left[K_z(z^{\prime }) - \mathcal {T}(z^{\prime }) - \mathcal {A}(z^{\prime }) \right] {\rm d}z^{\prime }. \end{equation}

(5)

However, this would mean that the derived value for σ_z(z) at some finite z would depend the properties of our mass model, tilt and axial terms not just in the range [0, z] but also for the range [0, ∞).

Equation (3) is numerically appealing to solve since, unlike many previous methods (e.g. Garbari, Read & Lake 2011; Garbari et al. 2012), it does not require any numerical differentiation and is therefore rather robust to noise in the data. Furthermore, equation (3) is valid for any gravity theory and can therefore be used as a constraint on alternative gravity models. In this sense, we follow the early pioneering work of Hill (1960) who attempted to also measure K_z directly without reference to the Poisson equation.

To connect the vertical acceleration K_z to the surface mass density of the disc, however, we must specify a gravitational model. For standard Newtonian weak field gravity, this is given by the Poisson equation:

\begin{equation} \nabla ^2 \Phi = \frac{\mathrm{\partial} ^2 \Phi }{\mathrm{\partial} z^2} + \hspace{-8.53581pt}\underbrace{\frac{1}{R}\frac{\mathrm{\partial} V_{\rm c}^2(R)}{\mathrm{\partial} R}}_{{\rm `rotation \,curve {'}\, term}:\,\,\mathcal {R}} \hspace{-8.53581pt}= 4\pi G\rho \end{equation}

(6)

where V_c(R) is the circular speed (rotation) curve at radius R, and ρ is the total matter density. Notice that for a flat rotation curve, V_c = const., and the ‘rotation curve’ term |$\mathcal {R}$| vanishes. If |$\mathcal {R}$| does not vanish, then it appears as a shift (potentially as a function of height z) to the recovered ρ_DM that can be corrected for at the end of the calculation (Garbari et al. 2012). Equation (6) can be rewritten as

\begin{equation} \frac{\mathrm{\partial} ^2 \Phi }{\mathrm{\partial} z^2} = 4 \pi G \rho (z)_{\rm eff} \end{equation}

(7)

where

\begin{equation} \rho (z)_{\rm eff} = \rho (z) - \frac{1}{4 \pi G R} {\frac{\mathrm{\partial} V_{\rm c}^2(R)}{\mathrm{\partial} R}}. \end{equation}

(8)

The contribution of this term to the mass density profile can be quantified via the Oort constants A and B (Binney & Tremaine 2008):

\begin{equation} \frac{1}{4 \pi G R} \frac{\mathrm{\partial} V_{\rm c}^2(R)}{\mathrm{\partial} R} = \frac{B^2 - A^2}{2 \pi G}. \end{equation}

(9)

Note that |B| < |A|, meaning that the terms in equation (9) are negative. Thus the effective density ρ(z)_eff will be higher with the inclusion of a non-zero rotation curve term. If the rotation curve term is erroneously neglected it will be absorbed into the density profile, yielding an over-estimate of the DM density, assuming the baryonic contribution is well constrained. Taking the most precise values for A and B from Table 1 (Mignard 2000) yields a value for equation (9) of ∼0.1 GeV/c², which is roughly a third of the expected local DM density (e.g. Read 2014). For any accurate measure of the local DM density derived from real data we will need to incorporate this correction, but we leave this for future work.

Table 1.

Open in new tab

Values of Oort constants. We include the F giants even though the errors for them are substantially larger to show that, within current uncertainties, the Oort constants A and B do not depend on stellar type.

Source	Stellar	A	B
	type	(km s⁻¹kpc⁻¹)
Mignard (2000)	K-M giants	14.5 ± 1.0	−11.5 ± 1.0
Branham (2010)	F giants	14.85 ± 7.47	−10.85 ± 6.83
Branham (2011)	G giants	14.05 ± 3.28	−9.30 ± 2.87

Source	Stellar	A	B
	type	(km s⁻¹kpc⁻¹)
Mignard (2000)	K-M giants	14.5 ± 1.0	−11.5 ± 1.0
Branham (2010)	F giants	14.85 ± 7.47	−10.85 ± 6.83
Branham (2011)	G giants	14.05 ± 3.28	−9.30 ± 2.87

Table 1.

Open in new tab

Values of Oort constants. We include the F giants even though the errors for them are substantially larger to show that, within current uncertainties, the Oort constants A and B do not depend on stellar type.

Source	Stellar	A	B
	type	(km s⁻¹kpc⁻¹)
Mignard (2000)	K-M giants	14.5 ± 1.0	−11.5 ± 1.0
Branham (2010)	F giants	14.85 ± 7.47	−10.85 ± 6.83
Branham (2011)	G giants	14.05 ± 3.28	−9.30 ± 2.87

Source	Stellar	A	B
	type	(km s⁻¹kpc⁻¹)
Mignard (2000)	K-M giants	14.5 ± 1.0	−11.5 ± 1.0
Branham (2010)	F giants	14.85 ± 7.47	−10.85 ± 6.83
Branham (2011)	G giants	14.05 ± 3.28	−9.30 ± 2.87

Neglecting |$\mathcal {R}$| and integrating both sides of equation (6), we derive the familiar result:

\begin{equation} \Sigma _z(z) = \frac{|K_z|}{2\pi G} \end{equation}

(10)

where Σ_z(z) is the surface mass density of the disc.

The overall flow of our method is now apparent. We model the tracer density ν(z), the mass density distribution ρ(z) = ρ_DM(z) + ρ_baryon(z) (which gives us the K_z term), the tilt term |$\mathcal {T}(z)$|⁠, and the axial term |$\mathcal {A}(z)$|⁠. These elements each have a number of parameters, and so in total each model will have an N-dimensional parameter space. Specifying values for each of these parameters will then give us quantitative values for each of these elements, which can then be inserted into equation (3) to derive |$\sigma _z^2(z)$| for that parameter space point. We then compare ν(z) and |$\sigma _z^2(z)$| [and σ_Rz(z) as part of the tilt term model] to data via a χ² test, and then change the values of our parameters.

Note that the only assumptions that go into equation (3) are that the motions of stars obey the collisionless Boltzmann equation, and that the galaxy is in dynamic equilibrium. The assumption of dynamic equilibrium can be negated somewhat by the presence of disequilibria (‘wobbling’) in the disc caused by e.g. the Sagittarius merger (Purcell et al. 2011; Gómez et al. 2013) or in part by the presence of spiral arms (Faure, Siebert & Famaey 2014). However, the impact of these effects on the measurement of the local DM density is estimated to be approximately 10 per cent (Widrow et al. 2012; Read 2014), less than the corrections arising from the tilt and rotation curve terms.

In practice, further assumptions are necessary to model the individual components on the RHS of equation (3). However, this method can accommodate almost any model for each of these elements – the only strict requirement is that each element can be defined at the z-values corresponding to the centres of the bins used to calculate the tracer density and velocity dispersion. Hence we call this method ‘non-parametric’ – we can in principle have many more parameters for each model than there are data points. In the following subsections we describe how we model each element for this particular study. In Section 2.6 we discuss model selection using the Bayesian evidence, which could potentially allow us to compare alternate models.

As the galaxy is close to being axisymmetric in both the thin disc (Dehnen 1998; Hogg et al. 2005; Aumer & Binney 2009; Bovy, Hogg & Roweis 2009; Pasetto et al. 2012b; Aghajani & Lindegren 2013) and the thick disc (Pasetto et al. 2012a; Sharma et al. 2014), the axial term |$\mathcal {A}(z)$| is expected to be small. For this study we assume complete axisymmetry, and thus take |$\mathcal {A}(z) = 0$|⁠. If the data show significant non-axisymmetry in the Milky Way, the axial term could be modelled in a similar way to the tilt term as described below.

Note that if the axial term were non-negligible, in the context of the Milky Way, this would imply the presence of significant spiral arms or a bar. Such features are time dependent, implying that if |$\mathcal {A}(z) \ne 0$|⁠, we should also worry about the time dependence of the gravitational potential Φ(t) and, therefore, also the partial time derivative of the distribution function ∂f/∂t (e.g. Monari, Famaey & Siebert 2016). It is not clear if there can exist a regime in which |$\mathcal {A}(z) \ne 0$| and such time dependence can be ignored. Here, we simply note that since empirically |$\mathcal {A}(z) \ll \mathcal {T}(z) \, \forall z$| (e.g. Pasetto et al. 2012a,b), to leading order we can also safely neglect the time dependence of the distribution function, as we have already done.

2.2 Tracer density model

To apply the Jeans equations we bin data to obtain the tracer density ν and the velocity dispersions σ_z and σ_Rz. Thus at a bare minimum we only have to define ν and σ_Rz, and derive σ_z via equation (3), at the bin centres, i.e. at a discrete set of z values.

For this work we model the tracer density as a sum of N exponential discs:

\begin{equation} \nu (z) = \sum _i \nu (0)_i \exp \left(-\frac{z}{h_i} \right), \end{equation}

(11)

where for the ith disc ν(0)_i is the tracer density at z = 0 and h_i is the scale height. The number of exponential discs can be increased until the Bayesian evidence shows no improvement (see Section 2.6). This method avoids overfitting as each disc is smooth, while still giving freedom to describe more complex data. From this point onwards we consider only z ≥ 0.

2.3 Baryon parametrization

For this study we use a simple two-parameter model to describe the baryon distribution:¹

\begin{equation} \rho _{\rm baryon}(z) = \frac{1}{4 \pi G} \left| \frac{K_{\rm bn}D_{\rm bn}^2}{(D_{\rm bn}^2 + z^2)^{1.5}} \right|. \end{equation}

(12)

The K_bn parameter sets the mass of the disc, and has dimensions of acceleration, while D_bn controls the scale height of the disc, and has dimensions of length. Expressed in terms of the K_z parameter of equation (3), the baryon profile becomes

\begin{equation} K_{z,{\rm baryon}} = - \left[ \frac{K_{\rm bn} z}{\sqrt{z^2 + D_{\rm bn}^2}} \right]. \end{equation}

(13)

While this model is not likely realistic for the Milky Way (Flynn et al. 2006; McKee, Parravano & Hollenbach 2015), it has been applied to observational data by Kuijken & Gilmore (1989a,b), and also more recently by Zhang et al. (2013). When applying our method to real data we will consider more sophisticated baryon models, but for this initial study where we are primarily interested in testing our methodology, equation (13) will suffice.

2.4 Dark matter parametrization

The simplest way to parametrize the DM density profile ρ_DM is to assume it is constant with z, as done in previous work (e.g. Bahcall 1984; Kuijken & Gilmore 1989a,b; Creze et al. 1998; Garbari et al. 2011, 2012). This assumption works well at low z: for a spherical NFW halo with a scale radius of 20 kpc the mid-plane value is correct within 10 per cent up to a height of z ∼ 3 kpc. For some analyses we also make this assumption, and set ρ_DM(z) = ρ_{DM, const.}.

However, this assumption does not allow for the exploration of several interesting effects such as a flattened, oblate halo or a DD. To accommodate such phenomena we add a DD on top of the constant DM. The DD is described using the same disc model as we use for the baryons, yielding a total DM profile:

\begin{equation} \rho _{\rm DM}(z) = \rho _{\rm DM, const.} + \frac{1}{4 \pi G} \left| \frac{K_{\rm DD}D_{\rm DD}^2}{(D_{\rm DD}^2 + z^2)^{1.5}} \right|. \end{equation}

(14)

To ensure this disc does not simply become degenerate with the baryonic disc we give the scale height D_DD a higher prior range than that for D_bn. In terms of the K_z parameter of equation (3), the DM profile is

\begin{equation} K_{z,{\rm DM}} = - \left[ 2Fz + \frac{K_{\rm DD} z}{\sqrt{z^2 + D_{\rm DD}^2}} \right] \end{equation}

(15)

where F = 2πGρ_{DM, const}. This disc parametrization of the non-constant DM profile can describe both the accreted DD and a flattened DM halo, and so henceforth we refer only to a ‘dark disc’.

Additional DD terms could in principle be added to ρ_DM(z) with parameters K_{DD, n} and D_{DD, n}, giving a total density profile of

\begin{equation} \rho _{\rm DM}(z) = \rho _{\rm DM, const.} + \frac{1}{4 \pi G} \sum _n \left| \frac{K_{{\rm DD},n}D_{{\rm DD},n}^2}{(D_{{\rm DD},n}^2 + z^2)^{1.5}} \right|. \end{equation}

(16)

The prior ranges on the nth DD parameters K_{DD, n} and D_{DD, n} can be set in relation to those of the (n − 1)^th DD, e.g. requiring K_{DD, n} < K_{DD, n − 1}, meaning each DD is less massive than the previous one. The number of DD to add could be determined from the data via the Bayesian evidence, with DD terms being added up until this degrades significantly. For this study, however, we limit ourselves to one DD term.

2.5 Tilt term

The tilt term from equation (2) links radial and vertical motion. The importance of this term has been noted previously, e.g. Kuijken & Gilmore (1989a), Smith, Wyn Evans & An (2009), Büdenbender et al. (2015). Here, we demonstrate that it is possible to deal with this term while remaining within our vertical, one-dimensional framework. Given the data quality currently available for the Milky Way we are required to make a well-motivated assumption about the radial form of |$\mathcal {T}$|⁠, but in the Gaia-era we will be able to directly measure this from the data.

We first assume that the radial profiles of the tracer density and the (R, z)-velocity dispersion are exponentials with scale radii of R₀ and R₁, respectively:

\begin{eqnarray} \nu (R,z) = \nu (z) \exp (-R/R_0), \end{eqnarray}

(17)

\begin{eqnarray} \sigma _{Rz} (R,z) = \sigma _{Rz}(z) \exp (-R/R_1). \end{eqnarray}

(18)

With this assumption the tilt term becomes

\begin{equation} \mathcal {T}(R,z) = \sigma _{Rz} (R,z) \left[\frac{1}{R} - \frac{1}{R_0} -\frac{1}{R_1} \right]. \end{equation}

(19)

If the disc were observed not to be an exponential then an alternative model could be easily applied at this stage. Indeed, similar but more complex models have featured previously in the literature (see e.g. Binney et al. 2014).

We then model the (R, z)-velocity dispersion as a power law at a given galactocentric radius R:

\begin{equation} \sigma _{Rz}(R,z) = \left. A\left(\frac{z}{{\rm kpc}}\right)^n \right|_R. \end{equation}

(20)

We take R = R_⊙, and also R₀ = R₁, simplifying the tilt term to a model described by the parameters A, n, and R₀:

\begin{equation} \mathcal {T}(\mathrm{R}_{\odot }, z) = \left. A \left(\frac{z}{{\rm kpc}}\right)^n \right|_{\mathrm{R}_{\odot }} \left[\frac{1}{\mathrm{R}_{\odot }} - \frac{2}{R_0} \right]. \end{equation}

(21)

Note that we are not affected by the assumption R₀ = R₁ since these two parameters are trivially degenerate in equation (21). It remains the case that any observational constraints on R₀ and/or R₁ can be used as priors on equation (21), where these would constrain the term 2/R₀. The description of σ_{R, z}(z) of equation (20) fulfils σ_{R, z}(z = 0) = 0 by construction and fits remarkably well to different populations, as shown in Fig. 1. This figure presents (z, σ_Rz) data points for high- and low-metallicity populations (blue and red points, respectively) from the SDSS/SEGUE survey as analysed and presented in Büdenbender et al. (2015), with a sign correction applied.² The high metallicity sample has |$\overline{[{\rm Fe/H}]} = -0.07$| and |$\overline{[\alpha {\rm /H}]} = 0.11$|⁠, while the low-metallicity sample has |$\overline{[{\rm Fe/H}]} = -0.89$| and |$\overline{[\alpha {\rm /H}]} = 0.34$|⁠. Metallicity can be used as a proxy for age, with the high-metallicity sample being younger than the older, low-metallicity population (e.g. Loebman et al. 2011).

$σRz for two populations of G-type dwarf stars in the solar neighbourhood from the SDSS/SEGUE survey (Büdenbender, van de Ven & Watkins 2015). The blue data points are from a younger, high-metallicity population, with $\overline{[{\rm Fe/H}]} = -0.07$ and $\overline{[\alpha {\rm /H}]} = 0.11$, while the data points in red are from an older, low-metallicity population with $\overline{[{\rm Fe/H}]} = -0.89$ and $\overline{[\alpha {\rm /H}]} = 0.34$. The blue and red lines are power laws in the form of equation (20) fitted to the blue and red data points, respectively. The older, lower metallicity stars probe further above the disc plane. Dashed lines indicate extrapolation from data.$

Figure 1.

σ_Rz for two populations of G-type dwarf stars in the solar neighbourhood from the SDSS/SEGUE survey (Büdenbender, van de Ven & Watkins 2015). The blue data points are from a younger, high-metallicity population, with |$\overline{[{\rm Fe/H}]} = -0.07$| and |$\overline{[\alpha {\rm /H}]} = 0.11$|⁠, while the data points in red are from an older, low-metallicity population with |$\overline{[{\rm Fe/H}]} = -0.89$| and |$\overline{[\alpha {\rm /H}]} = 0.34$|⁠. The blue and red lines are power laws in the form of equation (20) fitted to the blue and red data points, respectively. The older, lower metallicity stars probe further above the disc plane. Dashed lines indicate extrapolation from data.

Open in new tab Download slide

Taking these data points we fit power-law models as per equation (20), with parameters A = 123.99 km² s⁻², n = 1.16 for the high-metallicity population (blue) and A = 180.08 km² s⁻², n = 1.44 for the low-metallicity population (red). The low metallicity population samples further above the disc plane and populates the canonical thick disc of the Milky Way. Populations such as this are more interesting for local DM searches as they allow us to probe higher z regions where the baryon mass contribution begins to drop away leaving behind the DM halo. However as we go higher in the disc the tilt term becomes increasingly important, as illustrated by Fig. 2, which shows the variable ζ(z):

\begin{equation} \zeta (z) = \frac{K_{z,{\rm DM}}}{ K_{z,{\rm DM}} - \mathcal {T} } \end{equation}

(22)

where K_{z, DM} = −2Fz, the constant DM density term, with F = 267.65 km² kpc⁻² s⁻² and R₀ = R₁ = 2.5 kpc, the same values as are used to later generate our mock data in Section 3. Compared to the thin disc population (blue line), the effects of the tilt term become important at much lower z values. In short, to probe local DM we want to use thick disc stars probing higher-z ranges, but the cost is that we must include the tilt term in our analyses.

Figure 2.

The importance of tilt, as quantified by equation (22). The blue and red lines correspond to tilt terms derived from the high- and low-metallicity population fits of Fig. 1. The deviation arising from the tilt term increases more rapidly with height for the low-metallicity population (red line), which probes the high-z region most useful for determining the DM profile. Dashed lines indicate extrapolation from data. For this case R₀ = R₁ = 2.5 kpc, and F = 267.65 km² kpc⁻² s⁻².

Open in new tab Download slide

Ignoring the tilt term will always cause an underestimation of the local DM density, if all other components of the model such as baryons are held steady. Looking at equation (21), we note that R₀ < R_⊙ and A > 0, meaning the tilt term |$\mathcal {T}(\mathrm{R}_{\odot }, z)$| is always negative. Then considering equation (3) we see that to fit to |$\sigma _z^2(z)$| in the absence of the tilt term |$\mathcal {T}(z^{\prime })$|⁠, the K_z(z′) term, a negative term, is forced to become less negative in order to compensate. This requires a lower mass density, which, if the baryon density profile is held constant, manifests itself as a decrease in the DM mass density.

2.6 Statistical analysis and MultiNest

The model outlined above gives us an N-dimensional parameter space. To explore this parameter space and derive limits on the various observables we adopt nested sampling (Skilling 2006) as implemented in the publicly available MultiNest code (Feroz & Hobson 2008; Feroz et al. 2009; Feroz et al. 2013). MultiNest is a tool for Bayesian inference and parameter estimation. Bayes theorem is

\begin{equation} P(\theta | D, \mathcal {M}) = \frac{P(D | \theta , \mathcal {M}) P(\theta | \mathcal {M})}{P(D|\mathcal {M})} \end{equation}

(23)

where |$\mathcal {M}$| is the given model, θ is the set of parameters for that model, and D is the data. The left hand side, |$P(\theta | D, \mathcal {M})$|⁠, is known as the posterior, while the three terms on the right are the likelihood|$P(D | \theta , \mathcal {M}) = \mathcal {L}(\theta )$|⁠, the prior|$P(\theta | \mathcal {M})$|⁠, and the Bayesian evidence|$P(D|\mathcal {M})$|⁠.

The Bayesian evidence, a.k.a. the marginal or model likelihood, is a measure of how well the model performs given the data, and can be expressed as

\begin{equation} \mathcal {Z} = P(D|\mathcal {M}) = \int P(D|\theta ,\mathcal {M})P(\theta |\mathcal {M}) \,{\rm d}\theta . \end{equation}

(24)

The performance of two different models given the same data can be compared using the Bayes factor:

\begin{equation} B_{01} = \frac{P(D|\mathcal {M}_0)}{P(D|\mathcal {M}_1)}. \end{equation}

(25)

Assuming Gaussian errors it is possible to derive an empirical scale relating the Bayes factor to the strength of evidence for one model over another, as done in Trotta (2008). There, a |ln B₀₁| value of less than 1 is considered inconclusive, while values of 1.0, 2.5 and 5.0 are considered to give weak, moderate and strong evidence, respectively.

MultiNest takes a given prior probability distribution and likelihood function and calculates the posterior distribution and Bayesian evidence. Our likelihood function is based on the χ² distribution:

\begin{equation} \mathcal {L}(\theta ) = \exp {\left( \frac{-\chi ^2}{2} \right)}. \end{equation}

(26)

The value of χ² is simply

\begin{equation} \chi ^2 = \chi ^2_\nu + \chi ^2_{\sigma ^2_{z}} + \chi ^2_{\sigma _{Rz}}, \end{equation}

(27)

where

\begin{eqnarray} \chi ^2_\nu = \sum _j \frac{(\nu _{{\rm data},j} - \nu _{{\rm model},j})^2}{\mathit {SD}_{\nu ,j}^2}, \end{eqnarray}

(28)

\begin{eqnarray} \chi ^2_{\sigma _z} = \sum _j \frac{(\sigma _{z,{\rm data},j} -\sigma _{z,{\rm model},j})^2}{\mathit {SD}_{{\sigma ^2_z},j}^2}, \end{eqnarray}

(29)

\begin{eqnarray} \chi ^2_{\sigma _{Rz}} =\sum _j \frac{(\sigma _{Rz,{\rm data},j} -\sigma _{Rz,{\rm model},j})^2}{\mathit {SD}_{{\sigma _{Rz}},j}^2}. \end{eqnarray}

(30)

Note that if the reconstruction model does not contain a tilt term (e.g. |$\mathcal {T} = 0$|⁠) then |$\chi ^2_{\sigma _{Rz}}$| is set to zero.

Table 2 shows the prior ranges used for our analyses. We derive credible regions (CRs) for the DM density parameters by taking its posterior distribution and marginalizing over the remaining parameters. As outlined above our model has several components that can be turned on or off, such as the DD. Using the Bayesian evidence as calculated by MultiNest it is potentially possible to perform model comparison to determine which reconstruction model best fits the data. This idea will be explored in greater depth in subsequent studies.

Table 2.

Open in new tab

Prior ranges for parameters. Gaussian priors are described by a median μ and a standard deviation |$\mathit {SD}$|⁠. Note that ν(0) and σ_z(0) are the tracer density and velocity dispersion at z = 0, while quantities subscripted with 0, such as ν₀ and |$\mathit {SD}_{\nu , 0}$|⁠, are the values derived from data in the 0th bin, whose bin centre z₀ > 0. Tracer density ν(0) has units of [stars kpc⁻³]. K_bn and K_DD terms have units [km² kpc⁻¹ s⁻²].

Model parameter		Range or Gaussian μ & SD	Type
Baryons	K_bn	μ = 1500, \|$\mathit {SD}=150$\|	Gaussian
	D_bn	μ = 0.18 kpc, \|$\mathit {SD}=0.02$\| kpc	Gaussian
Constant DM	ρ_DM	[2, 20] × 10⁻³ M_⊙ pc⁻³	Linear
		[0.076, 0.796] GeV cm⁻³
DD	K_DD	[0, 1500]	Linear
	D_DD	[1.5, 3.5] kpc	Linear
Tilt term	A	[0, 200] km² s⁻¹	Linear
	n	[1.0, 1.9]	Linear
	R₀	μ = 2.5 kpc, \|$\mathit {SD}=0.5$\| kpc	Gaussian
Tracer density	ν(0)	[0, ν₀ + \|$2\times \mathit {SD}_{\nu , 0}$\|]	Linear
	h	[0.4, 1.4] kpc	Linear
Velocity disp.	σ_z(0)	[\|$\sigma _{z,0} - 2 \times \mathit {SD}_{\sigma _z, 0}$\|⁠,	Linear
		\|$\sigma _{z,0} + 2\times \mathit {SD}_{\sigma _z, 0}$\|] km s⁻¹

Model parameter		Range or Gaussian μ & SD	Type
Baryons	K_bn	μ = 1500, \|$\mathit {SD}=150$\|	Gaussian
	D_bn	μ = 0.18 kpc, \|$\mathit {SD}=0.02$\| kpc	Gaussian
Constant DM	ρ_DM	[2, 20] × 10⁻³ M_⊙ pc⁻³	Linear
		[0.076, 0.796] GeV cm⁻³
DD	K_DD	[0, 1500]	Linear
	D_DD	[1.5, 3.5] kpc	Linear
Tilt term	A	[0, 200] km² s⁻¹	Linear
	n	[1.0, 1.9]	Linear
	R₀	μ = 2.5 kpc, \|$\mathit {SD}=0.5$\| kpc	Gaussian
Tracer density	ν(0)	[0, ν₀ + \|$2\times \mathit {SD}_{\nu , 0}$\|]	Linear
	h	[0.4, 1.4] kpc	Linear
Velocity disp.	σ_z(0)	[\|$\sigma _{z,0} - 2 \times \mathit {SD}_{\sigma _z, 0}$\|⁠,	Linear
		\|$\sigma _{z,0} + 2\times \mathit {SD}_{\sigma _z, 0}$\|] km s⁻¹

Table 2.

Open in new tab

Prior ranges for parameters. Gaussian priors are described by a median μ and a standard deviation |$\mathit {SD}$|⁠. Note that ν(0) and σ_z(0) are the tracer density and velocity dispersion at z = 0, while quantities subscripted with 0, such as ν₀ and |$\mathit {SD}_{\nu , 0}$|⁠, are the values derived from data in the 0th bin, whose bin centre z₀ > 0. Tracer density ν(0) has units of [stars kpc⁻³]. K_bn and K_DD terms have units [km² kpc⁻¹ s⁻²].

Model parameter		Range or Gaussian μ & SD	Type
Baryons	K_bn	μ = 1500, \|$\mathit {SD}=150$\|	Gaussian
	D_bn	μ = 0.18 kpc, \|$\mathit {SD}=0.02$\| kpc	Gaussian
Constant DM	ρ_DM	[2, 20] × 10⁻³ M_⊙ pc⁻³	Linear
		[0.076, 0.796] GeV cm⁻³
DD	K_DD	[0, 1500]	Linear
	D_DD	[1.5, 3.5] kpc	Linear
Tilt term	A	[0, 200] km² s⁻¹	Linear
	n	[1.0, 1.9]	Linear
	R₀	μ = 2.5 kpc, \|$\mathit {SD}=0.5$\| kpc	Gaussian
Tracer density	ν(0)	[0, ν₀ + \|$2\times \mathit {SD}_{\nu , 0}$\|]	Linear
	h	[0.4, 1.4] kpc	Linear
Velocity disp.	σ_z(0)	[\|$\sigma _{z,0} - 2 \times \mathit {SD}_{\sigma _z, 0}$\|⁠,	Linear
		\|$\sigma _{z,0} + 2\times \mathit {SD}_{\sigma _z, 0}$\|] km s⁻¹

Model parameter		Range or Gaussian μ & SD	Type
Baryons	K_bn	μ = 1500, \|$\mathit {SD}=150$\|	Gaussian
	D_bn	μ = 0.18 kpc, \|$\mathit {SD}=0.02$\| kpc	Gaussian
Constant DM	ρ_DM	[2, 20] × 10⁻³ M_⊙ pc⁻³	Linear
		[0.076, 0.796] GeV cm⁻³
DD	K_DD	[0, 1500]	Linear
	D_DD	[1.5, 3.5] kpc	Linear
Tilt term	A	[0, 200] km² s⁻¹	Linear
	n	[1.0, 1.9]	Linear
	R₀	μ = 2.5 kpc, \|$\mathit {SD}=0.5$\| kpc	Gaussian
Tracer density	ν(0)	[0, ν₀ + \|$2\times \mathit {SD}_{\nu , 0}$\|]	Linear
	h	[0.4, 1.4] kpc	Linear
Velocity disp.	σ_z(0)	[\|$\sigma _{z,0} - 2 \times \mathit {SD}_{\sigma _z, 0}$\|⁠,	Linear
		\|$\sigma _{z,0} + 2\times \mathit {SD}_{\sigma _z, 0}$\|] km s⁻¹

3 MOCK DATA SETS

In this paper we apply our method only to mock data in order to hone and verify it. This mock data is ‘as good as it gets’, in the sense that it has no measurement errors, nor observational biases added to it, and is drawn from relatively simple, known distribution functions. While dynamically realistic ‘N-body’ mocks are preferred, it has already been shown that 1D Jeans analyses are robust to the presence of local non-axisymmetric structure in the disc (Garbari et al. 2011). Furthermore, N-body mocks are expensive; even state-of-the-art simulations do not approach the local sampling expected from Gaia. Finally, the expense of building well-sampled N-body mocks makes it challenging to explore a large parameter space of models, including models with and without tilt, or with and without a DD. For these reasons, we focus here on simpler mock data, similarly to the Read (2014) review. We will consider more dynamically realistic mocks, and mocks that give a faithful representation of the expected Gaia data uncertainties in future work.

We generate here a variety of mock data sets as outlined in Table 3, with different samplings (10⁴, 10⁵, and 10⁶ tracer stars), with and without a tilt term, and also with either no DD, a DD (dd), or a more massive ‘big dark disc’ (bdd). We assume the axial term |$\mathcal {A}$| and the rotation curve term |$\mathcal {R}$| are zero. Our simple baryonic disc model is set up to mimic the real Milky Way, with a scale height parameter of D_bn = 0.18 kpc and a surface density of Σ_bn = 55.53 M_⊙ pc⁻², similar to those measured near the solar position (Flynn et al. 2006; Read 2014; McKee et al. 2015). The value of the F parameter (see equation 15) corresponds to a DM density of ρ_{DM, const} = 10 × 10⁻³ M_⊙ pc⁻³ = 0.38 GeV cm⁻³. For each scenario we generate 20 mock data sets, allowing us to explore the effects of Poisson noise over a range of realizations.

Table 3.

Open in new tab

Mock data parameters. _X is the number of stars sampled, e.g. 10⁴, 10⁵, 10⁶, and _M is the mock number, ranging from 0 to 19. Empty spaces indicate that a certain mock data set does not include that element. The baryon model corresponds to a baryonic surface density of Σ_bn = 55.53 M_⊙ pc⁻², while the F parameter corresponds to a constant DM density of ρ_{DM, const} = 10 × 10⁻³ M_⊙ pc⁻³ = 0.38 GeV cm⁻³. K_bn and K_DD terms have units [km² kpc⁻¹ s⁻²], while F has units [km² kpc⁻² s⁻²].

The mock data consist of a list of stars each with three pieces of data: the position z, the vertical velocity v_z, and the product of the vertical and radial velocities, v_Rv_Z. The first element, z, is generated by drawing randomly from an exponential tracer distribution with scale height h:

\begin{equation} \nu (z) = \exp \left(-\frac{z}{h} \right). \end{equation}

(31)

This gives us our list of stellar positions.

To derive the velocities for the mock catalogue we first must define a mass model and tilt. Taking the same parametrizations as described in Sections 2.3, 2.4, and 2.5 for baryons, DM, and tilt, respectively, we set their parameters as per Table 3. This allows us to calculate K_z and |$\mathcal {T}(z)$|⁠. Using equation (5) and its associated assumptions we can calculate C. We then use equation (3) to derive σ_z(z). For each star we take its position z′, find the value of σ_z(z′), and then draw a velocity from a Gaussian centred on v_z = 0 and with variance |$\sigma _z^2(z^{\prime })$|⁠.

To generate v_Rv_Z mock data, when necessary, we take the A, n, and R₀ parameters listed in Table 3 and calculate the σ_Rz(z) profile via equation (20). Taking each star's position z′, we calculate σ_Rz(z′), and draw a value of v_Rv_z from a Gaussian centred on v_Rv_z = 0 and with variance of |$\sigma _{Rz}^2(z^{\prime })$|⁠.

4 RESULTS

Here we present the results of our scans. We first investigate how the precision of the reconstruction varies with different numbers of tracer stars. We then look at the effects of the tilt term and the DD. For certain parts of these analyses the model used to generate the mock data and the model used to reconstruct the mock data are not the same – this is done to investigate the effects of ignoring terms (as in the case of the tilt term) or using incorrect models.

In the following figures we plot the marginalized posterior for the DM density ρ_DM(z), showing 68 per cent, 95 per cent, and 99.7 per cent CRs in dark, medium and light shading, respectively. The mock data profile is shown by the solid black line, while the median of the posterior distribution is shown by a solid line in the same colour as the CR. Binning of stars to calculate ν(z), σ_z(z), and σ_Rz(z) is performed such that each bin contains an equal number of stars. For this study we use 20 bins, and in Appendix A we briefly explore the effects of changing the number of bins used. For plots with constant DM density in both mock data set and reconstruction model we plot all 20 mock data sets together, while for non-constant DM density in either mock or reconstruction we show one representative figure in this section, and show the full set of figures in Appendix B (Figs B1 to B6).

Our method and code is set up to describe the tracer density as a sum of exponentials. To determine the necessary number of exponentials required to describe the data we can use the Bayes factor as described earlier in Section 2.6. Test reconstructions give a Bayes factor of 1.5 when comparing one exponential to two, and 3.3 when comparing one exponential to three, with one exponential favoured in both cases. This is to be expected given that the mock data are generated using a single exponential. Given this result we use a single exponential for all subsequent reconstructions.

4.1 Sampling

Fig. 3 shows reconstructions of the local DM density for using varying numbers of tracer stars. The mock data sets used here are the most simple case, thick_X_M as described in Table 3, and has no DD or tilt added. As expected the CRs shrink as the number of tracer stars is increased from 10⁴ stars up to 10⁶ stars. The SDSS has sampling of the order of 10⁴ stars (Zhang et al. 2013; Büdenbender et al. 2015) while data from Gaia will give upwards of 10⁶ tracer stars.

Figure 3.

Determination of the DM density profile for varying numbers of tracer stars. These plots show marginalized posteriors for ρ_DM(z) = ρ_{DM, const} for the 20 mock data sets generated for each sampling level. Dark, medium, and light shading show the 68 per cent, 95 per cent, and 99.7 per cent CRs, respectively. Green, purple, and orange colouring indicates that the 68 per cent, 95 per cent, or 99.7 per cent CR, respectively, contains the correct answer. The median value of each posterior is shown by a solid line in green, purple, or orange, while the DM density value used to generate the mock data is shown as a solid black line across the entire plot. The mock data and reconstruction models used contain no tilt term or DD terms. As the number of tracer stars is increased from 10⁴ to 10⁶ the CRs for ρ_{DM, const} naturally shrink around the mock data value.

Open in new tab Download slide

4.2 Tilt

In Fig. 4 we explore the effects of the tilt term. The left-hand set of CRs in Fig. 4 is the same as the right hand set of CRs from Fig. 3, with no tilt term in the mock data or reconstruction. In the centre set of Fig. 4, however, the mock data thick_tilt_1e6_M has the tilt term turned on, but the analyses are performed with the tilt term set to zero. This illustrates the danger of ignoring tilt as discussed earlier in Section 2.5: the reconstructions return narrow CRs, but as expected they systematically underestimate the value of ρ_DM, with the true ρ_DM lying outside even the 99.7 per cent CRs. This underestimation, however, is remedied by turning on the tilt term in the model, as shown in the right hand set CRs in Fig. 4, where the correct DM density is always at least within the 95 per cent CR. The inclusion of extra parameters to describe the tilt term does, however, increase the size of the CRs.

Figure 4.

Determination of the DM density profile using 10⁶ tracer stars and exploring the effects of including or neglecting the tilt term in the mock data sets and reconstruction models. These plots show marginalized posteriors for ρ_DM(z) = ρ_{DM, const} for the 20 mock data sets generated for each tilt scenario. Dark, medium, and light shading show the 68 per cent, 95 per cent, and 99.7 per cent CRs, respectively. Green, purple, and orange colouring indicates that the 68 per cent, 95 per cent, or 99.7 per cent CR, respectively, contains the correct answer, while pink colouring indicated that the correct answer lies outside even the 99.7 per cent CR. The median value of each posterior is shown by a solid line in green, purple, or orange, while the DM density value used to generate the mock data is shown as a solid black line across the entire plot. The left set of mocks shows the same result as the right set of Fig. 3, with no tilt term in mock data or reconstruction model. The centre set has a tilt term in the mock data but not in the reconstruction model, yielding a systemic underestimation of ρ_DM. The right set has a tilt term in both mock data and reconstruction model, demonstrating that our method can successfully deal with the tilt term.

Open in new tab Download slide

Just as our determination of ρ_DM(z) is dependent on the tilt term, the tilt term is in turn dependent on its input parameters, A, n, and R₀ = R₁. While we have been able to fit for A and n using G-type dwarf data (Section 2.5), for R₀ we have taken the canonical value of R₀ = 2.5 ± 0.5 kpc from Binney & Tremaine (2008), and further made the assumption that σ_Rz(R, z) has the same radial profile as ν(R, z) (i.e. R₀ = R₁). When using only a single population, determination of R₀ and R₁ will be important, as illustrated in Fig. 5. This figure shows the 2D marginalized posterior for ρ_{DM, const} and R₀ for a reconstruction of mock data set thick_tilt_1e6_0 with a model including tilt, and demonstrates the degeneracy that exists between ρ_{DM, const} and R₀ = R₁. If using multiple tracer population, each will have different R₀ and/or R₁, but will all have their motions dictated by the same potential. This will help us break the degeneracy between R₀, R₁, and ρ_DM. Further, with the advent of Gaia data we will be able to directly measure the radial profile of σ_Rz and ν(R, z) for a given set of tracer stars.

2D marginalized posterior of the constant DM density ρDM(z) = ρDM, const and the R0 parameter of the tilt term (equation 19), generated from mock data set thick_tilt_1e6_0 reconstructed using a model with tilt. Contours show the 68 per cent, 95 per cent, and 99.7 per cent CRs. Marginalization and plotting performed using Barrett (Sebastian Liem, private communication).

Figure 5.

2D marginalized posterior of the constant DM density ρ_DM(z) = ρ_{DM, const} and the R₀ parameter of the tilt term (equation 19), generated from mock data set thick_tilt_1e6_0 reconstructed using a model with tilt. Contours show the 68 per cent, 95 per cent, and 99.7 per cent CRs. Marginalization and plotting performed using Barrett (Sebastian Liem, private communication).

Open in new tab Download slide

4.3 Dark disc

Fig. 6 shows mock data set thick_1e6_0-19, thick_dd_1e6_0, and thick_bdd_1e6_0 reconstructed using models with and without a DD component. The top left panel shows the same CRs as seen in the left hand set of Fig. 4 and the right hand set of Fig. 3.

Figure 6.

Determination of the DM density profile using 10⁶ tracer stars and exploring the effects of including or neglecting a DD in the mock data sets and reconstruction models. For comparison the top left panel shows the same reconstructions as the right hand set from Fig. 3, i.e. basic mock data sets with no DD, reconstructed with a constant DM density. The remaining panels each show one representative mock and reconstruction, with the full set of mocks and reconstructions given in Appendix B. These panels show marginalized posteriors for ρ_DM(z), with dark, medium, and light shading indicating the 68 per cent, 95 per cent, and 99.7 per cent CRs, respectively. The median value of the posterior is shown by the solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line. From left to right the columns show reconstructions of mock data sets containing no DD (ρ_{DM, const} only), a DD, or a ‘big’ dark disc (BDD). The top centre and top right panels show the determination of a constant DM density from the DD and BDD mock data sets, respectively. The bottom row shows the reconstruction of each of the mock data sets using a model with a constant DM term and a DD term. Light grey vertical lines indicate the bin centres.

Open in new tab Download slide

The left column of Fig. 6 shows the reconstruction of mock data sets with a constant DM density profile: mocks 0–19 in top left and mock 0 in bottom left. The reconstruction in the top panel uses a model with a constant DM density, while the bottom panel uses a model with an additional DD component. The DD reconstruction exhibits a disc structure even though the correct answer is constant ρ_DM. This is likely due to a bias in the hyper-volume set by the priors – the prior range on the DD parameters goes between no DD (K_DD = 0) and a maximal DD (K_DD = 1500), and thus the bulk of the parameter space features at least some DD. There is no ‘negative dark disc’ to counteract this effect and push the mean of the prior range back to no DD.

In the centre and right columns of Fig. 6 we reconstruct mock data sets with a DD of different masses: thick_dd_1e6_0 and thick_bdd_1e6_0 (the ‘big dark disc’). A constant DM density reconstruction (top row) is able to contain the thick_dd_1e6 DD within the 95 per cent CR almost to the last bin, but fails to contain the big dark disc beyond z = 1.3 kpc. Adding a DD term allows the reconstruction to fit to the mock data DM profile nicely, as shown in the bottom centre and bottom right panels of Fig. 6.

When working with real data, however, we will not have the luxury of knowing the DM profile underlying the data. The black mock data model line will not be there. In subsequent studies we will explore the potential use of the Bayesian evidence, as calculated by MultiNest, to determine if a DD is justified by observational data.

4.4 Tilt and dark disc

Here, we combine the two elements discussed in previous sections, the tilt term and the DD. Fig. 7 shows reconstructions of the mock data set thick_bdd_tilt_1e6_0. In the top panel the reconstruction model contains neither DD nor tilt term. Again, we see the same effects as we did previously. The missing tilt term yields a consistent underestimation of the DM density, and the constant DM density envelope is too narrow to encompass the density range of the DD. The consistent underestimation can be remedied by adding a tilt term to the reconstruction model, as shown in the middle panel, where the reconstruction model has a tilt term included. Both problems can be resolved in tandem by using a reconstruction model with both tilt and DD terms, as shown in the bottom panel of Fig. 7.

Figure 7.

Determination of the DM density profile using 10⁶ tracer stars with a combination of a tilt term and a DD in the mock data (thick_bdd_tilt_1e6_0) and using a variety of reconstruction models. These plots show marginalized posteriors for ρ_DM(z), with dark, medium, and light shading indicating the 68 per cent, 95 per cent, and 99.7 per cent CRs, respectively. The median value of the posterior is shown as the solid blue line, while the DM density profile used to generate the mock data is shown as a solid black line.

Open in new tab Download slide

5 CONCLUSIONS

In this paper we have presented a new method of determining the vertical DM density profile, and thus the local DM density. The key equation of this method, equation (3), depends only on the assumption of dynamical equilibrium. In practice, further assumptions are made to describe the components going in to equation (3); however the scope of possible models is wide, and in principle model selection using the Bayesian evidence can be used to determine the best model. We leave exploration of this last aspect to future work.

The importance of the tilt term has been previously noted (Kuijken & Gilmore 1989a; Smith et al. 2009; Büdenbender et al. 2015), and to derive an accurate value for the local DM density it is vital that our method be able to deal with this term. The baryonic contribution to the galactic density profile drops rapidly as z increases, leaving DM as the dominant component. Thus to probe the DM density profile sampling thick disc stars that travel higher above the plane are preferable. For these populations the tilt term is more important, as illustrated in Fig. 2. Failure to include the tilt term in the analysis leads to a systematic underestimation of the local DM density, as explained in Section 2.5 and demonstrated in Section 4.2.

One of the novel aspects of our method is that it can deal with the tilt term while remaining within the confines of the one-dimensional z-direction Jeans equation, which can be seen in Fig. 4. With only the data currently available for the Milky Way, this requires several well-motivated assumptions, as described in Section 2.5. However, with data from Gaia we will be able to directly measure the radial profile of the tilt and tracer density, removing the need for such assumptions. While for this paper we have disregarded the rotation curve term (equation 6), we note that an accurate determination of this will be necessary for the implementation of this or other z-direction methods to real data.

Non-spherical DM density profiles, such as oblate haloes or accreted DDs, can also be fitted using our method by incorporating a DD term into the reconstruction model, which is shown in Fig. 6. Our method can also reconstruct mock data sets containing both a tilt term and a DD, as shown in Fig. 7.

We thank Michael Feyereisen, Sebastian Liem, Pat Scott, Roberto Trotta, Stephen Feeney, Roberto Ruiz de Austri, Scott Tremaine, Vera Gluscevic, and Juna Kollmeier for very useful discussions. SS acknowledges support from the Swedish Research Council (VR, Contract No. 350-2012-268). JIR would like to acknowledge support from STFC consolidated grant ST/M000990/1 and the MERAC foundation. GB acknowledges support from the European Research Council through the ERC starting grant WIMPs Kairos.

1

In the context of the values of the prior ranges given in Table 2, the gravitational constant equals |$G = 4.229\times 10^{-6}\, {\rm km}^2\, {\rm kpc}\, {\rm M}_{\odot }^{-1}\, {\rm s}^{-2}$|⁠.

2

Private communication with authors.

REFERENCES

Aghajani

T.

Lindegren

L.

2013

A&A

551

A9

Month:	Total Views:
December 2016	2
March 2017	6
April 2017	6
May 2017	3
June 2017	4
August 2017	7
September 2017	2
November 2017	2
December 2017	2
January 2018	9
February 2018	3
March 2018	4
April 2018	3
May 2018	26
July 2018	3
August 2018	9
September 2018	3
November 2018	8
January 2019	1
March 2019	1
April 2019	8
May 2019	10
June 2019	2
July 2019	12
August 2019	7
September 2019	9
October 2019	4
November 2019	10
December 2019	9
January 2020	11
February 2020	7
March 2020	6
April 2020	6
May 2020	2
June 2020	4
July 2020	2
August 2020	5
September 2020	18
October 2020	6
December 2020	6

Article Contents

A non-parametric method for measuring the local dark matter density Free

Abstract

1 INTRODUCTION

2 METHOD

2.1 Deriving a general 1D jeans method

2.2 Tracer density model

2.3 Baryon parametrization

2.4 Dark matter parametrization

2.5 Tilt term

2.6 Statistical analysis and MultiNest

3 MOCK DATA SETS

4 RESULTS

4.1 Sampling

4.2 Tilt

4.3 Dark disc

4.4 Tilt and dark disc

5 CONCLUSIONS

REFERENCES

APPENDIX A: VARIATION DUE TO THE NUMBER OF BINS

APPENDIX B: Dark Disc Reconstructions with Multiple Mocks

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only

A non-parametric method for measuring the local dark matter density